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SUMMARY 


The steady form of the full potential equation, in conservation form, is employed to 
analyze and design a wide variety of complex aerodynamic shapes. The nonlinear method 
is based on the theory of characteristic signal propagation coupled with novel flux biasing 
concepts and body-fitted mapping procedures. The resulting code is vectorized for the 
CRAY-XMP and the VPS-32 supercomputers. 

Use of the full potential nonlinear theory is demonstrated for a single-point supersonic 
wing design and a multipoint design for transonic maneuver/supersonic cruise/maneuver 
conditions. Achievement of high aerodynamic efficiency through numerical design is veri- 
fied by wind tunnel tests. Other studies reported here include analyses of a canard/wiag/ 
nacelle fighter geometry. 
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1.0 INTRODUCTION 


Development of supersonic/hypersonic configurations has traditionally relied on 
linearized methods and hypersonic impact theory. These approaches can treat complex 
geometries with minimum response time and cost, providing wide data coverage in terms of 
Mach number, angle of attack, trim deflection, yaw angle, etc. Shortcomings are present, 
however, in both the impact and the linearized methods. For the former, interference be- 
tween surface elements is totally ignored in implementations such as classical Newtonian, 
tangent wedge and cone theories. Crossflow interactions and stagnation point singularities 
are also implicitly disregarded. In the latter, shocks, vorticity, and entropy wakes and 
layers are excluded. Furthermore, superposition of elementary solutions such as those for 
thickness and angle of attack freely used in linear models are, strictly speaking, invalid at 
supersonic/hypersonic speeds. 

Modern vehicle concepts such as the Advanced Tactical Fighter (ATF) attempt to 
achieve an effective compromise between the transonic maneuver and supersonic cruise/ 
maneuver conditions. Multiple design considerations of this type impose stringent con- 
straints on the aerodynamic shape of the vehicle to achieve high buffet-free lift perfor- 
mance with reduced trim drag. The analysis and design of these modern vehicle shapes 
is a tremendous challenge, requiring increasingly sophisticated nonlinear methods ranging 
from the full potential theory to the Navier-Stokes equations. 

Full potential approximations include enough physics of the flow to allow realistic 
optimization and permit consideration of mutual interference of highly integrated, closely 
coupled arrangements to provide improvements to aerodynamic efficiencies achievable using 
linear methodology. This approach, in conjunction with the use of modern high-speed 
computers, achieves the objective of economic computational design that is responsive to 
conceptual aircraft development efforts. 

This report presents a state-of-the-art technique to solve the steady form of the full 
potential equation employing several novel concepts such as 1) flux biasing for capturing 
shocks, 2) implicit approximate factorization scheme for computational efficiency, 3) wake 
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treatment, 4) numerical mapping for treatment of complex geometries, and 5) vector cod- 
ing for supercomputers. The steady form of the full potential equation is mainly used for 
treating predominantly supersonic flows with embedded subsonic regions, while the un- 
steady full potential equation is employed to treat blunt-nosed configurations, transonic 
Mach number flows and time-accurate unsteady simulations (oscillating wings). A unified 
full potential method 1 is currently under development which can handle flows across the 
Mach number range (subsonic, transonic, supersonic), including the supersonic marching 
technique 2-6 as a special case. 

This report demonstrates the use of this full potential nonlinear aerodynamic pre- 
diction methodology to treat a wide variety of complex configurations, including a single 
point and a multipoint design of advanced fighter wings. Results shown include geometries 
with canard, nacelle and vertical tail. 

In designing a fighter wing, linear theory 7,8 is used to establish candidate optimum 
thickness and zero drag-due-to-lift twist, camber and variable camber deflections at super- 
sonic speeds. Nonlinear potential flow analyses are employed to capture embedded shock 
waves at transonic 9,10 and supersonic conditions 2-6 and then subsequently weaken the 
wave system through parametric redesign. Boundary layer analysis 11 follows the inviscid 
design/analysis to assess the flow quality of the nonlinear potential design. The extent 
of trailing edge separation in particular is evaluated. The general approach is schemati- 
cally indicated in Fig. 1 and represents a summary of the numerical design experience at 
Rockwell covering the HiMAT, forward swept wing, SAAB and Air Force/Navy Research 
Technology contract studies. 

The full potential code employed in this paper is operational on several computer sys- 
tems, such as the CYBER 176 and VAX serial machines and the CRAY-XMP and VPS-32 
supercomputers. Analysis of a complete fighter-like configuration requires 500 seconds on 
the CYBER 176 and about 20-35 seconds on the VPS-32 or CRAY-XMP class machines. 
The execution time required to run a case depends on the number of marching plane 
calculations and on the number of grid points at each marching plane. A typical fighter 
calculation may involve 400 marching planes to cover the entire length of the configuration 
with an average of 75 x 20 grid points per plane. 
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A brief description of the full potential methodology is presented here for completeness. 
More details can be found in Refs. I through 6 and 9. 


INITIAL DESIGN. DESIGN CRITIQUE/REDESIGN, DESIGN CRITIQUE. 

3-D LINEAR 3-D NONLINEAR 3-D BOUNDARY LAYER 


CANDIDATE 


CANDIDATE 



•DIRECT DESIGN -ACCURATE ANALYSIS •PROBLEM DETECTION 

•QUICK TURNAROUND -PROBLEM DETECTION ‘BOUNDARY LAYER 

•RELATIVELY INEXPENSIVE ‘SHOCK OCCURRENCES SEPARATION 

• CRITICAL PRESSURES 


Fig. 1. Numerical design approach. 
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2.0 METHODOLOGY 


As mentioned earlier, the steady form of the full potential equation is used for treating 
predominantly supersonic flows with pockets of subsonic regions. 


2.1 Treatment of the Steady Equation 


The steady, conservative full potential equation cast in an arbitrary coordinate system 
defined by f x,y,z ), rj = r](x,y,z) and £ = ((x,y,z) can be written as 




( 1 ) 


where the density p is given by 
P = 

The nature of Eq. (1) can be analyzed by studying the eigenvalue system of Eq. (1) 
combined with the irrotationality condition in the and (£,£) planes. A detailed 

discussion on this and the nomenclature can be found in Ref. 4. Therefore, only the final 
results are presented here. 

1. At a grid point, the marching direction ( is hyperbolic (an — (U 2 fa 2 )) < 0 and the 
total velocity q is supersonic, q > a. This point will use the algorithm of Ref. 3. The 
quantity a u is (f 2 + f 2 + f 2 ). 

2. At a grid point, the marching direction f is elliptic (an — (U 2 /a 2 )) > 0, but the total 
velocity q is supersonic, q > a. This point will be treated by a transonic operator 
with a built-in density biasing based on the magnitude of [l — (a 2 /g 2 )]. This case is 
termed Marching Subsonic Region (MSR). 

3. At a grid point, the direction f is elliptic (an — [U 2 fa 2 )) > 0 and the total velocity 
q is subsonic, q < a. This point will be treated by a subsonic central-differenced 
operator. This case is called Total Subsonic Region (TSR). 


-(¥) 


1 1 /( 7 ~ 1 ) 


MioiU+'+Vh+Wh-l) 


( 2 ) 
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Figure 2 shows the schematic of a fuselage-canopy forebody geometry with an em- 
bedded MSR and TSR present in the supersonic flow. To solve this problem, the marching 
scheme of Ref. 3 is used when [an — (f/ 2 /a 2 )] is positive. First, march from the nose 
up to the plane denoted by (A-B) in Fig. 2 using the method of Ref. 3. Then, between 
(A-B) and (C-D), which bound the subsonic bubble (MSR and TSR), use a relaxation 
scheme and iterate until the subsonic bubble is fully captured. Then, resume the marching 
scheme from the plane (C-D) downstream of the body. For blunt-nosed configurations, 
the unsteady method of Ref. 9 is used to generate the starting blunt body solution. 


Treatment of d/d$ (py) Term 


At a grid point (t + l, j,k), the derivative in the marching direction f is given by 


d 




i 

supersonic 


+ <*-*«>£(> j)^ 

marching subsonic 


( 3 ) 
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Fig. 2. Embedded subsonic bubble in a supersonic flow. 
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where 


d refers to backward differencing 
d refers to forward differencing 

&i = 1 if [an — (U 2 / a 2 )] <0 (supersonic with respect to f) 

= 0 if [on - (tf 2 /a 2 )] > 0 (subsonic with respect to f) 
and i + 1 is the current marching plane. 

The first term in Eq. (3) corresponds to the supersonic marching operator and the 
second term is the subsonic operator. By using a local linearization procedure, Eq. (3) can 
be expressed in terms of <p only. The supersonic operator is given by 


d_ 





+ 


( UV \ d aj. ( UW \ d K± TT 

" 15- J y - + (“« - -r j + u ‘ 


A<f> = - <f>i) 


( 4 ) 


and the subsonic relaxation operator is given by 


± (~U\ j, d_ 

3s yj) i+l 3s 


«? 


(on <f>f +CIl2^ij +Gi3^),- + 1 


( 5 ) 


where 

P?+t = Pi+i ~ * (Pi+i ~ P?) i for U > 0 
v = max 

The superscript (n + 1) denotes the current relaxation cycle for a subsonic bubble 
calculation. 
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Treatment of djdr\ (py) Term 


d_ 

dri 


«)- -4 W) 


y+i/2 


o a, 


( 6 ) 


3+1/2 


supersonic 


marching subsonic 


When 6{ + i = 1, that is, the point is supersonic with respect to f, only the first term 
in Eq. (6) is used and the biased density p is defined by (for V > 0) 


Pj+ 1/2 “ (i “ ^+1/2) Pj+ 1/2 + 2 U 3+ 1 / 2 + P*3~ 1 ) > 


( 7 ) 


where P = max [0, 1 - 022 (o 2 /V 2 )]. 

In Eq. (7), the evaluation of p* depends on whether the flow is conical or nonconical. 
For conical flows, all p* quantities are evaluated at the i th plane. For nonconical flows, 
at each nonconical marching plane, initially p* is set to be the value at the : th plane and 
then subsequently iterated to convergence by setting p* to the previous iterated value of 
p at the current t + 1 plane. 

When the point is elliptic in the marching direction, the density biasing p is based on 
”= max [0, 1 — (a 2 /g 2 )]. 

Combining the various terms of Eq. (1) as represented by Eqs. (3)— (7) together with 
the terms arising from [p(W/J)f] will result in a fully implicit model. This is solved 
using an approximate factorization implicit scheme 4 . Details on the initial and boundary 
conditions, wake treatment, and geometry and grid setup can be found in Ref. 5. 

For treatment of complex geometries, the body-fitted coordinate system (r, f,» 7 , f) is 
generated using the procedure outlined in Ref. 12. 

2.2 Treatment of Combined Yaw and Angle of Attack 

A complete analysis of an aircraft configuration must address the vehicle’s performance 
under angle of attack (a) and sideslip (/ 3 ) flight conditions. For asymmetric configurations, 
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a similar problem exists in that the flow field is not symmetric, and therefore requires a 
solution of the entire crossflow plane. 

Given a symmetric configuration at a sideslip angle of /? = 0, only the half-plane 
problem needs to be solved with the plane of symmetry boundary conditions imposed 
along K = 2 and (KM AX — 1), as shown in Fig. 3a (see Ref. 4 for a complete description 
of the nomenclature). Imposing the flow conditions along K = 1 to be the same a s the 
ones along K = 3 gives a tridiagonal system of equations for the L( operator that can 
be easily solved. For a symmetric configuration, when yaw (or sideslip) angle is present, 
the entire crossflow plane needs to be solved, as shown in Fig. 3b. In this case, the flow 
conditions along K = 1 are set to be the same as the ones along K = (KM AX — 2). This 
destroys the tridiagonal nature of the L ^ operator. A special routine has been developed 
to invert a matrix of the following type: 

rX X X 01 


XXX — 0 0 

XXX 0 



Lo x x x\ 


In the current formulation, positive angle of attack a represents a positive Cartesian 
velocity v in the freestream, and similarly positive yaw /? produces a positive w in the 
freestream. When both angle of attack and yaw are present, first the freestream is turned 


by an angle /? and then by a. Let ( x , y, z) be the inertial Cartesian system. After an initial 
yaw turn /?, let the wind axis system be (x', y\z'), and after a subsequent a turn let it 


become (z,y,z). Now, referring to Fig. 3c, 
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cos a 

sin a 

o' 


cos/9 0 

sin /? 
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— sin a 
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0 1 
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cos a cos / 3 

sin a 

cos a sin /? 


X 



r= 

— sin a 

cos/? 

cos a 

— sin a sin/? 


y 





-sin/? 

0 


cos/? 


z 



( 9 ) 


The freestream is now along x. The normalized freestream velocity potential is given by 


4L = — 


x cos a cos /? + y sin a + z cos a sin /?. 


( 10 ) 
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K-1 2 3 



•) PLANE OF SYMMETRY 
(YAW ANGLE * 0) 

f|.1 * *i.3 

♦j.KMAX - Pj. KM AX- 2 


KMAX-1 



b) PERIODIC CONDITION FOR 
YAW TREATMENT 

<j.1 * Pj.KMAX-2 
<j. KMAX * Pj.3 


Fig. 3a, b. Boundary conditions imposed for cases with and without yaw treatment. 


Y.v y.v 





Fig. 3c. Coordinate transformation for yaw and anglc-of-attack treatment. 
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Using Eq. (9), the drag, lift and side forces are easily calculated: 


Drag = F x cos a cos /? + F y sin a + F z cos a sin /? 

Lift = — F x sin a cos/? + F y cos a — F z sin a sin/? . (11) 

Side Force = -F x sin /? + F z cos / 3 

In Appendix B, the demonstration of the code for a combined yaw and angle of attack over a 
fighter-like configuration is presented to familiarize the user with the various input/output 
aspects of the code. 

More results on the yaw capability of the code can be found in Ref. 13. 

2.3 Nonisentropic Flows 

The full potential equation is based on the assumption that the flow is isentropic. 
This is true only for very weak shocks (Mach number normal to the shock is less than 1.4). 
Some researchers have tried to modify the full potential formulation to correct the results 
for entropy effects 14 . No such attempt is made in the present work. Instead, an Euler 
capability is being developed 15 (included in Appendix C) to treat flows with strong shocks 
and rotational effects. A sample result from this Euler development is provided in the 
Results section of this report. 

2.4 Salient Features of the Marching Code 

The marching code developed under this NASA contract is named Supersonic Implicit 
Marching Program (SIMP) and is available from COSMIC under the designation LAR- 
13413 16 . The code structure is described in Appendix A. Some of the salient features of 
the code are: 

• Equation in conservation form 

• Flux linearized upwind differencing in the marching direction 

• Conservative switch operators to treat embedded subsonic zones 
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• Treatment of wakes 


• Yaw and angle of attack 

• Complex geometry treatment (fuselage, canopy, wing, canard, nacelle, tail, multibody, 
etc.) 

• Numerical grid generation with constraints 

• Use of GEMPAK 17 or CDS 18 to generate geometry input files 

• Vectorized code for supercomputers. 
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3.0 APPLICATIONS 


To demonstrate the capability of this nonlinear full potential methodology, the appli- 
cation to several configurations is discussed: 

1. Nonlinear aerodynamic wing design for an advanced fighter configuration. 

2. Canard/wing/nacelle supersonic fighter analysis. " 

3. Shuttle, Transatmospheric Vehicle (TAV), and supersonic inlet cowl design. 

Also, in Appendix B — User’s Guide, a complete test case is presented for an advanced 
fighter configuration for a combined yaw and angle of attack flight condition. 

3.1 Case 1 — Nonlinear Aerodynamic Wing Design 

A cooperative effort between Rockwell International and the National Aeronautics and 
Space Administration-Langley Research Center was conducted to determine the effect on 
supersonic aerodynamic characteristics of increasing wing sweep on a North American 
Rockwell fighter design 19,20 . The effort was also to provide validation of the nonlinear full 
potential analysis code described in this report. 

The configuration used in the study is a preliminary design version of a Rockwell 
fighter concept (see Fig. 4). This concept includes a slender forebody and prominent canopy 
blended smoothly into a highly swept wing center section. The outboard wing panel, 
selected for transonic performance requirements, extends outward from about 37 percent 
of the semispan and has a 48° swept leading edge. Leading- and trailing-edge devices 
for high lift and roll control extend along most of the span of the outer wing panel. The 
wing is twisted and cambered for transonic maneuvering. The propulsion system consists 
of two engines located beneath the center-wing section in nacelles which are blended into 
the lower surface inboard of the 37-percent semispan. The nozzles are vectorable in pitch. 
Twin vertical tails are located at the outboard edge of the center-wing section and are 
canted out 20°. Control surfaces located at the trailing edge of the wing center section, in 
conjunction with the canted vertical tails and vectorable nozzles, provide pitch control. 
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Fig. 4. General arrangement. 



A model of the configuration was tested in the Langley Unitary Plan Wind Tunnel 
and Rockwell’s Trisonic Test Facility. In addition to the 48° leading-edge outboard wing 
panel of the concept, two 55° leading-edge panels, and two new 48° leading-edge panels 
were designed and tested. One of the 55° panels was twisted and cambered for a Mach 1.6 
maneuver design point. The two new 48° panels were designed with leading- and trailing- 
edge flap systems to best meet the subsonic/transonic/supersonic cruise and meneuver 
design points. For details of the wing panel design procedure, see Ref. 20. Aerodynamic 
force data for the 48° baseline panel and the 55° panels can be found in Ref. 19. Experi- 
mental data (aerodynamic force and surface pressure) for the new 48° panels and surface 
pressure data for the 55° panel will be reported in a forthcoming NASA report. 

The computational model of the wing-body-tail-nacelle fighter under development is 
shown in terms of surface grid plots in Figs. 5 and 6. A typical marching plane contained 
a 75 x 20 mesh generated by the elliptic grid solver of Ref. 12. Figure 7 shows the cross- 
sectional grid and circumferential surface pressures at various axial stations in front of the 
nacelle face at the cruise flight conditions. Figures 8 and 9 show results at two different 
axial stations where the nacelles are present. The presence of a shock around the nacelles 
is clearly seen in Fig. 8. Details of the crossflow velocity vectors (projection of the total 
velocity vector on a unit sphere whose center is at the apex of the fighter configuration) 
around the wing-body-nacelle geometry is shown in Fig. 9. In Fig. 9, the freestream is 
aligned predominantly along the x-axis and away from the geometry and is seen as inflow 
crossflow velocity vectors. Crossflow velocity information provides insight into the behavior 
of vortical singularities and plays a key role in activating the density biasing switches in 
the code necessary to simulate crossflow shocks. 

Linearized theory and full potential estimates of the 55° twisted and cambered wing 
panel configuration are presented in Fig. 10. Fully turbulent skin friction drag for a 
mean aerodynamic chord test Reynolds number of 1.56 x 10® is used for this assessment. 
Examination of the results indicates the design is an aerodynamically efficient one, taking 
into consideration nominal scale effects. 

Figure 11 shows comparisons of chordwise pressures at two different span stations 
(60% and 80% span). The results show that the full potential predictions are in good 
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Fig. 5. Surface grid for ATF with nacelle. 




Fig. 7. Pressure distribution on a fighter-like (ATF) configuration; M , » = 1.6, a = 1.24. 


17 




PRESSURE CONTOURS GRID 


Fig. 8. Pressure contours and grid for ATF with nacelle; M w = 1.6, a = 5°, f = 0.46. 



Fig. 9. Pressure contours and grid for ATF with nacelle; = 1.6, a = 5°, f = 0.56. 
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Turbulent Skin Friction 
Cd f = 0.0129 unitary , = 1.56 x 10° 

Linear Analysis 

a = 4.46 deg 
C L = 0.32 

C Dp = C Dw + C* Dt - 0.109 + 0.02185 = 0.03275 
C M = -0.061 

L/D = 0.32/(0.0129 + 0.03275) = 7.0 

Full Potential 

a = 4.46 deg 
C L =0.311 
C Dp = 0.0325 
C M = -0.0579 

L/D = 0.311/(0.0129 + 0.0325) = 6.85 
“Cambered plate fuselage 

Fig. 10. Pretest M = 1.6 maneuver point design drag assessment. 
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Fig. 11. 




; Moo = 1.6, a = 1.24. 











agreement with Rockwell’s experimental data. Figure 12 shows comparison of overall forces 
and moments in terms of Cl, Cd p and Cm- 

Supersonic cruise and maneuver pretest assessment of the new 48° leading-edge design 
is summarized in Fig. 13. Linear predictions indicate lift-drag ratio levels of 3.25 and 6.19, 
respectively, for the proposed Unitary Plan Wind Tunnel test condition. Comparison with 
the 55° point design result of Fig. 10 indicates a 13% reduction results from multipoint 
compromises associated with wing sweep and airfoil leading edge radius. Two-thirds of 
the penalty is associated with thickness considerations and one-third with drag due to lift. 
Figure 14 shows the lift, drag and pitching moment results for the 48° leading-edge sweep 
multipoint design. The full potential predictions are in good agreement. 

The impact of multipoint considerations on the nacelle off untrimmed lift-drag ratio 
is presented in Fig. 15. For the 48° wing configuration, a 7.6% reduction at the nominal 
design point results from the decrease in planform sweep and increase in airfoil leading edge 
radius, twist and camber. This penalty is modest and overstated by full potential analysis 
which is slightly optimistic for the subsonic edge (55°) case and somewhat pessimistic for 
the supersonic edge (48°). 

3.2 Case 2 — Canard- Wing-Nacelle Fighter Analysis 

The full potential method of Refs. 2-6 can handle extremely complex geometries. This 
is demonstrated by applying the method to analyze the complex geometry 21 of Fig. 16, 
which has a canard, clipped wing tip, canopy and a flow-through nacelle mounted on the 
undersurface of the fuselage. Figure 17 shows the computational geometry and surface 
gridding. Note that the boundary layer diverter and swept nacelle were not modelled in 
the computations. Computations were performed for this configuration at Moo = 2 and 
a = 4°. Figures 18 and 19 show the crossflow streamlines, surface pressures, pressure 
contours and crossflow velocity vectors at two different model axial stations. Figure 18 
shows results at an axial station where the fuselage, wing, canard wake and canard are 
all present. The nodal singularity in pressure contour that is present at lower wing-body 
junction regions corresponds to a saddle singularity of crossflow streamlines, as shown in 
Fig. 18. Note the pressure match along the canard wake cut, Fig. 18. Figure 19 shows 
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Fig. 12. Comparison of measurement with predictions at M = 1.6 for nonlinear point 
design, A le =55 deg, nacelle off. 
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Turbulent Skin Friction 
Cd f = 0.0129 unitary, Rn c = 1.56 x 10 fl 

Linear Analysis 

a = 1.24 deg 
C L = 0.1 

C Dp = C Dw + C* Dl = 0.0149 + 0.0031 C Dp 
= 0.180 
C M = -0.0136 
L/D = 3.24 

Full Potential 

a = 1.24 deg a = 5.22 

C L = 0.119 C L = 0.336 

C Dp = 0.0164 C Dp = 0.04075 

C M = -0.0284 C M = -0.051 

L/D = 4.06 L/D = 6.26 

‘Cambered plate fuselage 

Fig. 13. Pretest M = 1.6 multipoint design drag assessment. 


a = 5.22 
C L = 0.32 

= C Dw + C* Dl = 0.0149 + 0.0239 
= 0.0388 
C M - -0.035 
L/D = 6.19 
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Fig. 14. Comparison of measurement prediction at M = 1.6 for first multipoint design, 
A le = 48 deg, nacelle off. 
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Fig. 15. Impact of sweep on lift-drag ratio at M = 1.6, nacelle off. 
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Fig. 16. Langley canard-wing fighter configuration. 



Fig. 17. Computational geometry and surface gridding for Langley fighter configuration. 
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(I)- NODAL SINGULARITY 


D 


Fig. 18. Solution for Langley fighter configuration; = 2.0, a = 4.0, x = 11.0; 
(a) streamline, (b) surface pressure coefficient, (c) pressure contour, (d) veloc- 
ity contour. 
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Fig. 19. Solution for Langley fighter configuration; Moo — 2.0, a = 4.0, x = 14.0; 

(a) streamline, (b) surface pressure coefficient, (c) pressure contour, (d) veloc- 
ity contour. 
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results at a station where the nacelle is present. The formation of a shock around the 
nacelle is clearly seen in Fig. 19. The upper and lower center plane pressure contours at 
Moo = 2.0 and a = 4° are shown in Fig. 20. The bow shock, canopy shock, nacelle shock 
and expansion wave are evident in this figure. 

3.3 Case 3 — Miscellaneous Cases 

3.3.1 Shuttle Aerodynamics 

Appendix C contains the reprint of a paper 6 presented in January 1985 (AIAA-85- 
0272) that deals with applications of the full potential method to three dimensional geome- 
tries including multibody configurations. The paper presents the results of calculations for 
the isolated Shuttle Orbiter and the mated Shuttle configuration (Orbiter, External Tank, 
and Solid Rocket Boosters). 

3.3.2 Transatmospheric Vehicle 

The paper of Ref. 6 also presented some preliminary results for analysis of a pro- 
posed transatmospheric vehicle (TAV) concept. Figure 21 presents results of additional 
calculations on an alternate TAV configuration. Analysis of these TAV concepts indicates 
that the rapid execution times of the vectorized full potential code makes it a very useful 
design/analysis tool for preliminary design of TAV configurations. 

Figures 22-24 show application of the full potential code to a TAV-like configuration 
with fuselage-mounted vertical tails. 

3.3.3 Supersonic Inlet Cowl Design 

Appendix C includes a paper 15 given in July 1985 (AIAA-85-1703) which illustrates 
the use of the SIMP program for the design of a twisted-cone inlet spike. The objective 
was to determine the proper location of the swept cowl lip such that the oblique shock 
emanating from the compression cone is attached at the cowl lip at the design supersonic 
Mach number. The SIMP calculations were compared with the results of an Euler solver 15 
under development and, as expected, the SIMP code overpredicts the pressure on the upper 
surface where a nonisentropic shock is formed. This problem illustrates the limitations of 
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the SIMP code and shows the need for an Euler method for supersonic flow computations 
which must capture strong shock waves. 



Fig. 20. Upper and lower centerplane pressure contour for Langley fighter configuration; 
Moo = 2.0, a = 4.0. 
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Fig. 23. Gridding 
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4.0 CONCLUSIONS 


This report demonstrates the use of a nonlinear full potential aerodynamic prediction 
capability developed at Rockwell International under a contract from NASA-Langiey Re- 
search Center. The method is now routinely employed to analyze and assist in the design 
of advanced fighter wings at supersonic flight conditions. Also, the nonlinear method is be- 
ing used to analyze extremely complex geometries that include complete fighter geometries 
(canard, nacelle, vertical tail and wake effects). Use of supercomputers and vector coding 
make application of this powerful nonlinear method very attractive and cost-effective. 

Appendices are included which give details of the code structure, input data, a sample 
test case and user’s guide, and relevant publications. 
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APPENDIX A — CODE STRUCTURE 


CODE ORGANIZATION 

The SIMP analysis code is applicable to arbitrary wing-body-nacelle-tail arrange- 
ments from moderate supersonic Mach numbers ( M , ~ 1.2) to values of the hypersonic 
similarity parameter MS < 1. The lower code limit is governed by the extent of the 
embedded subsonic flow while the upper limit results from a breakdown in the isentropic 
assumption for strong shock waves. Also, since the potential theory is irrotational, the 
modeling of any vortices is not attempted. 

The program is written in FORTRAN V language. It can be executed on any CDC 
machine (CYBER 176, CDC 7600), as well as on the CRAY-XMP and CYBER 205. For 
a cross-plane ( rj , £) grid of 25 x 140, the program requires 660,000 words of memory. The 
program consists of a main routine and several subroutines. A brief description of the code 
along with input instructions needed to execute the code are given in this Appendix. 

Program MAIN 

Program MAIN coordinates the entire operation. A flowchart describing the various 
operations performed by the MAIN program is given in Fig. Al. The MAIN program 
sets up the initial (known) data plane and the body-fitted grid system and performs the 
L{ and L n operators to advance the solution. The marching step size Af can either be 
prescribed or computed at each marching plane from a given Courant number and the 
maximum eigenvalue. The various read and write tapes used in the calculation are listed 
below. 

Program MAIN (Tape 1, Tape 2, Tape 3, Tape 4, Tape 5, Tape 7, Tape 8, Output, 

Tape 6 = Output). 

Tape 1: Output solutions for plot. 

Tape 2: Output solutions for restart. 
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Tape 3: Read in starting solutions for restart. 

Tape 4: Output solutions for restart. 

Tape 5: Input data. 

Tape 6: Solution output. 

Tape 7: Read tape containing solutions for subsonic region. 

Tape 8: Write tape for subsonic bubble calculation. 

Subroutine INVSXI and Subroutine INVSETA 

The factored implicit scheme for the governing full potential equation can be written 
as 

and it is implemented as follows: 

L ( {A*)*=R , L„(A<f>) = (A*)' 

+ &<(>• 

The subroutine INVSXI performs the f inversion, and the rj inversion is solved in subrou- 
tine INVSETA. 

Subroutine EIGEN (EIGENY, EIGENZ) 

This subroutine computes the maximum eigenvalue EIGENY in the (f , T]) plane and 
the maximum eigenvalue EIGENZ in the (f , £) plane. The expression used for calculating 
the eigenvalue is given in Ref. 5. The maximum eigenvalue information is then used to 
compute a marching step-size Af for a specified Courant number. 
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Subroutine NFORCE (PX, PY, PM, AREA, KG) 


At the end of each marching plane calculation, this subroutine computes the axial 
force, PX, vertical force, PY, and the side force, PZ, by integrating the pressure force acting 
on an elemental area, dA. The elemental area, dA, is computed from the transformation 
matrix using the formula (at a body point j = 2). 


dA = { [y { z ( - z f y^] 2 + [x(Z { - x f z*] 2 + [y*x f - y f x^] 2 } / d$d£. 
KG = 0, conical or blunt body nose force calculation 
= 1, rest of the body force calculation. 


The program also prints the force coefficients, Cl, Cp, and Cs (side force coefficient) 
information based on a prescribed reference area, and moment coefficients, Cm, Cy, and 
Cr about a given reference point ( Xq , Yo, Zq). Cy is the yawing moment coefficient 
and Cr is the rolling moment coefficient. Cs, Cy, and Cr are opposite in sign from the 
standard convention, since the W component of velocity is defined to be positive from the 
centerline out towards the left wing. 


Subroutine GEOM (NO, NRP) 


N9 = 0, 
> 0 , 


geometry data at X\ and are read in 

i 

geometry data at X\ is updated and X 2 is read in 


NRP = 0, constant z marching plane geometry calculation 
= 1, spherical marching plane geometry calculation 

Subroutine GEOM sets up the body grid points from a prescribed geometry shape. 
From the input geometry points, a key point system is established using cubic splines. 
These key points are then joined from one prescribed geometry station to the next to 
provide the geometry at any intermediate marching plane 12 . 
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Subroutine GRID 


Once the body points are obtained at a marching plane from GEOM, subroutine 
GRID sets up the entire crossflow plane grid using an elliptic grid solver that satisfies 
certain grid constraints. 

Subroutine METRIC 

This subroutine computes all the necessary transformation metrics and Jacobians at 
various node and half-node locations as required by the solution algorithm (L{ and L n 
operators). 

Subroutine UVW 

This subroutine computes all the contravariant velocities, U, V, and W, and the 
density p. 

Subroutine RHOBIAS 

This subroutine performs the density biasing in the (»/, £) plane based on characteristic 
signal propagation theory. This operation is essential to treat crossflow supersonic regions 
and to capture shock waves. 

INPUT DATA 

The input data can be divided into four parts: (1) header data describing mesh 
information, Mach number, angle of attack, aerodynamic coefficient reference quantities, 
center of gravity location, etc.; (2) detailed geometric coordinates defining configuration 
cross plane contours; (3) program update file directives defining code modifications; and 
(4) job control directives defining program and input/output file allocations (this depends 
on the particular computer). 
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Header Data 


A typical analysis of a complete configuration requires several regions of marching 
calculations for a complete analysis. Each region calculation has a different set of header 
instructions for describing grid parameters, wake information if pertinent, restart direc- 
tions, and number of mesh points for each patch of the region. A sample input is given 
in Appendix B for the configuration of Fig. A2, and a brief description of each variable is 
given in this section. 

Card* Symbol Format Description 

100 NMAX 15 Number of axial marching steps. 

If NMAX = 0, and ZTA1 = f and TAPER = F 
the code generates geometry and grid data 
at x = f for plotting. For NMAX^ 0, 
the code will march for NMAX steps unless 
XEND is encountered first. NMAX must 
include NCON iterations if applicable. 

(NMAX = 0 option for grid plot is provided 
to allow the user to review the quality 
of grid at various axial stations before 
the flow solver is turned on.) 

110 JMAX 15 Mesh points in the normal direction (17). 

Present maximum is 25. This can be 
increased by increasing the dimension. 

This number includes the J = 1 dummy line 
inside the geometry. 

120 KM AX 15 Mesh points in circumferential direction (£) 

(maximum value: 140). This number includes 

* Card number (e.g., Card 100) is not part of the input file. 
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the K = 1 and K =KMAX dummy lines. If this 

number is incorrectly specified, the code 

will reset KMAX properly using the 

KMAX = 1 + (IPT1-1) + (IPT2-1) + (IPT3-1) 

+ ■ • • + (IPTn-1) + 1 (definition of IPT 

follows in the next section), “n” is the 

number of patches. 

130 NRM 15 Number of grid regions (separated by 

dashed lines in Fig. A3). Maximum of 6 
allowed. 

140 NUO 15 Not used. 

150 NP 15 Selected surface output for every 

NP'steps. 

160 KWKMIN 15 K value of starting point of a patch 

containing wake (Fig. A3). 

170 KWKMAX 15 K value of ending point of a patch 

containing wake (Fig. A3). 

180 NCON 15 Number of iterations for conical starting 

solution (usually set to 30). To establish 
this starting solution, the geometry is 
initially assumed to be conical. The geometry 
at ZTAl is projected forward conically 
to a point at (0,0,0). (The nose of the 
configuration is assumed to be at (0,0,0). 

If the geometry is not input this way, 
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shift the geometry using PTNOSE and 
YSHIFT.) The solution is then obtained 
for this conical geometry based on 
NCON iterations. The conical solution 
is used as a starting solution for 
the nonconical case, beginning at 
ZTA1. The marching step size Af 
for the conical calculation is based 
on a specified Courant number (CFLIN). 

The user should be aware that NCON is 
included in the NMAX total. Also, ZTAl 
output values have no physical 
significance during conical calculation. 

190 NITER 15 Number of iterations to generate the 

marching grid using an elliptic grid solver. 
Usually set to 30. If grid routine fails, 
set this to 0 to analyze the geometry 
and the initial grid generated before grid 
relaxation (this is for debugging purposes). 
Set NITER back to 30 for flow field 
analysis. NMAX should be set to zero 
for analyzing the grid quality. 

200 NSPTI 15 Number of f locations for detailed 

flow field output (maximum of 10 locations). 

210 ITERGE 15 Number of global iterations for subsonic 

region calculations. 
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The marching calculation can encounter subsonic 
(MSR and TSR) regions, especially at low 
freestream supersonic Mach numbers ( M ~ 1.2 
to 1.8). For the first pass through subsonic 
regions, the relaxation routine in the code 
assumes sonic flux conditions (Ref. 4). If the 
extent of the subsonic region (usually confined 
near the body surface, especially near the 
leading edge) is very small, just the first 
pass with assumed sonic flux conditions will 
suffice for the marching calculation. However, 
if the subsonic zone is expected to be large, 
(around bumps on the geometry which might 
create detached shocks such as the canopy), 
several relaxation passes through the 
subsonic zone are essential for correct 
representation of the flow. ITERGE represents 
the number of relaxation passes for subsonic 
calculation. The user may not know ahead 
of time if a large subsonic zone is to be 
encountered during the marching calculation. 
The output prints the location of subsonic 
points encountered during the first pass. 

If too many subsonic points are predicted 
by the code during the first pass, then set 
ITERGE > 1 and TAPE8W = TRUE. 

Then, subsonic relaxation is carried out 
over the entire NMAX marching planes for 
ITERGE times. The relaxation solution from 
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each pass is stored on tape (TAPE8W) 
to aid in the relaxation process (Ref. 4). 

If the initial prescription of ITERGE 

value is not enough for satisfactory 

subsonic convergence (RMS change in 

density between two relaxation cycles 

should be a small value, preferably 10 -4 ), 

then additional subsonic iterations can 

be performed through the restart option 

(TAPER = TRUE, TAPE8W = TRUE) by reading 

the previously stored solution on TAPE7. 

During subsonic calculations, the marching 
step size A^ is to be kept constant. 

220 CFLIN F10.5 The CFL number. 

If DZTAIN (Card 230) is negative, the 
axial step size A^ is generated based 
on CFL number. The relationship between Af 
and CFL number through the maximum eigenvalue 
is given in Ref. 3. When the geometry change 
in the axial direction is minimal (nearly 
conical shape), the marching step size Af 
is given by the CFL number (usually CFLIN 
set to 5). If the geometry changes are very 
abrupt (emergence of wing, canopy, tail, or 
any other component) or drastically 
nonconical, then Af is prescribed by 
the user (see Fig. A4). 
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230 DZTAIN F10.5 Initial step size. For nonconical geometry 

calculations, DZTAIN is chosen to be either 
DZMIN or DZMAX. If DZTAIN is set to less 
than DZMAX, then during marching calculation, 
Af will be slowly increased to DZMAX. 

240 DZMAX F10.5 Maximum step size. 

250 DZMIN F10.5 Minimum step size. 

(DZMAX and DZMIN depend on the complexity 
of the geometry. Suggested value: 

DZMAX = total length/400 and 

DZMIN = DZMAX/2.) If DZMIN is set equal 

to DZMAX, then constant step size is used. 


260 

FSM 

F10.5 

270 

ALFA 

F10.5 

280 

THTO 

F10.5 


Freestream Mach number. 

Angle of attack (degrees). 

Angle of outer boundary (degrees). 

This angle must be larger than the bow 
shock wave in order for the code to capture 
the bow shock. Often the best way to 
choose this value is to calculate the 
bow shock wave half angle and add 10°. 
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290 


DZTA 


F10.5 First step size for the marching 

calculation after conical starting solution. 
Usually, this value is set to DZMIN. 

Af 


CONICAL 
PLANE 

300 ZTAl F10.5 Starting f location. If TAPER = 

TRUE, this value is overwritten by stored 
restart value. 

310 XEND F10.5 Final f location for this run. 

320 AMUl F10.5 1: first order accuracy in marching 

direction. 

0: second-order accuracy in marching 
direction. 

330 AMU2 F10.5 0: first-order accuracy in marching 

direction. 
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1: second-order accuracy in marching 
direction. 

Usually, first-order accuracy is used for rapidly 
varying geometries. For nearly conical cases, 
second-order accuracy is recommended. 

340 XWAKE F10.5 Wake starting location in the axial 

direction (see Fig. A3). 

350 BETANG F10.5 Angle of yaw (degrees). 

/ 

360 CHL F10.5 Geometry scale factor. If set to total 

length, f will be scaled from 0 to 1. 

If set to 1, actual dimensions of the 
geometry are used. Use of dimensional 
(CHL = 1) or nondimensional (CHL = i) 
option is left to user’s choice. 

370 PTNOSE F10.5 Axial geometry shift. Equal to negative 

of apex of the forebody (i.e., shifts 
configuration nose to ( = 0). 

x 

f=x+ PTNOSE 

380 YSHIFT F10.5 Vertical geometry shift (i.e., shifts 

configuration nose to tj = 0). 
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390 

XO 

F10.5 

Moment reference x location (unit~length) 

400 

YO 

F10.5 

Moment reference y location (unit~length) 

410 

AAA 

F10.5 

Reference area to compute aerodynamic 
force coefficients (unit~le;'gth 2 ). 

420 

ALL 

F10.5 

Reference length to compute aerodynamic 


moment coefficients (unit~length). 

XO, YO, AAA, and ALL are to be chosen 

(dimensional or nondimensional) based on CHL. 

430 OMEGA F10.5 Overrelaxation parameter for grid generation. 

Suggested value: 

1.0 (for vectorized code) 

1.75 (for scalar code). 

440 YAW L5 T: Calculation with yaw (full cross- 

plane grid). 

F: Without yaw (half-plane grid). 

450 NUGRID L5 T: Numerical grid generation (normally 

used). 

F: User must adapt code for his particular need. 

460 IRE AD L5 T: Read body geometry input which must 

be supplied in the format described in 
the next section titled “Geometry Data”. 
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F: Analytic geometry (which must be supplied 
by the user and inserted in subroutine 
GRID). 

470 RPLANE L5 T: Spherical plane marching (spherical 

plane marching is exercised only for 
conical flow calculations). 

F: Constant x plane marching. 

480 TAPER L5 T: Restart the calculation. 

F: Start the calculation from freestream 
<(> values. 

490 TAPEW L5 T: Write restart data on Tapes 2 and 4. 

F: No data storage for restart. 

500 TAPE8W L5 T: Write entire flow field data for 

subsonic iterations on Tape 8. 

F: No flow field data saved. 

When a tape read or tape write is set TRUE, 
the user will have to provide the necessary 
job control cards. 

510 FORCE L5 T: Compute aerodynamic forces and moments. 

F: No force computation. 

520 THTU(5) 515 Grid region terminal points ( k ) 

(see Fig. A3). These values are the K values 
of the points where the dashed lines intersect 
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the body. 


530 INU(5) 5F10.4 Polar angle (degrees) at respective 

terminal point. 

540 ZTAPT(IO) 10F8.4 x-locations of detailed flow field 

printouts. 

550 ISC 15 Number of patches (geometry) that define 

the cross-sectional shape of the configuration 
for this region of the configuration 
(see Figs. A3 and A5). (Maximum number 
of patches = 15.) 

560 NPT(15) 1515 Number of output points on each patch 

(maximum number of points per patch is 30). 

Geometry Data 

The cross-sectional geometry of a typical aircraft changes considerably in the axial 
direction due to emergence of various components such as canopy, wing, nacelle, and tail, 
etc. The marching computation, as it sweeps along the marching direction f, has to account 
for this geometry variation to set up the proper body-fitted coordinate system to aid in 
the application of body boundary conditions. To treat complex geometry cross sections, 
patches are introduced to define the geometry as indicated in Fig. A2. Using patches, a 
configuration is defined by several regions of cross sections. The number of patches defining 
a section is constant for a given region (Fig. A2). 

A complete computation over a configuration such as the one in Fig. A2 is usually done 
in segments rather than in one shot. The calculation starts from the nose and proceeds 
along Even within a region (defined by the same number of patches), the calculation 
might be done in segments using the restart option in the code. Restart is used any time 
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the calculation is halted and then continued with another run that picks up where the 
previous run left off. Pure restart is performed only when there is no alteration to the 
number of points along tj (JMAX) and along ( (KM AX), and no change in the number 
of grid points per patch between the previous run and the current restart run. If there is 
any alteration to the grid structure, the restart run will automatically perform a respace 
operation to interpolate the solution from the previous solution grid to the current grid. 
Respace is used whenever the following situations are encountered: 

1) Number of patches defining the cross section is changed. This situation occurs when 
the cross-sectional geometry becomes more complex. This is illustrated in Fig. A2. 

2) Number of JMAX and/or KMAX points is changed (even if the number of patches 
defining the cross section is kept the same as before). This situation often occurs for 
cases where a patch length is increasing with f. For example, a swept wing is very 
small when it first appears in the cross section of the geometry and only requires a 
few grid points for accurate computation of the flow field. However, as the analysis is 
continued in the $ direction, the wing patches grow' and will require more points for 
accurate flow field analysis. 

3) Number of grid points per patch is changed (even if KMAX is kept the same as before). 

Any time a respace is required, the code must be stopped. The code will automatically 
do a respace if KMAX or JMAX is different from the previous values of KMAX and JMAX. 

One may be able to compute the entire configuration using the same number of patches 
and same KMAX and JMAX values throughout to avoid the respace requirement. This will 
mean even in the forebody region of a configuration, where the cross-sectional geometry 
is usually simple, more grid points and more patches are to be used than necessary to 
adequately resolve the flow field. Use of the same number of patches and grid points 
for throughout the length of the configuration is generally not recommended. This can 
substantially increase the total execution time. 

Transitioning from one region to the next (number of patches is changed) requires an 
overlapped zone, as illustrated in Fig. A6, to allow for increased or decreased number of 
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patches in the next region. The extent of this overlapped zone must be sufficient to include 
at least the final three marching data planes of the prior region. In the overlapped region, 
the data from the previous region is interpolated onto the grids of the new region. For the 
example of Fig. A6, the results from the 3-patch region are interpolated onto a 4-patch 
region grid at the same x location. This is required in order to continue marching along 
the body with the new patch definition. 

Figure A6 illustrates how to transition from a fuselage computation to a wing-fuselage 
computation. First, the calculation is performed for the fuselage section denoted by 
REGION 1 which ends just prior to the starting point of the wing. This calculation might 
involve, say, three patches. Then, to introduce the wing, a four patch representation is 
used in REGION2. In the overlapped zone, the fuselage which is defined using a three 
patch representation in REGION 1 is represented by a four patch representation as part 
of REGION2. The second and third patch locations on the fuselage in REGION2 within 
the overlapped zone are chosen in the vicinity of where the leading edge of the wing is 
expected to emerge from the fuselage. 

Wake Geometry 

Behind the trailing edge of a lifting surface a cut is introduced (see Fig. A3), across 
which potential <f> jumps are imposed (the <f> jumps are computed at the trailing edge) 
to preserve density continuity across the flow through cut. Mathematical details of this 
so-called “wake model” are given in Ref. 4. The treatment of wake cut within the code 
requires the knowledge of starting and ending K index values of the upper wake cut and 
the lower one. Depending on the sweep of the trailing edge, the wake cut is appropriately 
modeled. This is illustrated in Fig. A3. The user has to define the shape of the trailing 
edge and also the starting x value where the wake begins to appear in the cross-sectional 
geometry (XWAKE). The wake cut is part of a patch which contains the wing also as 
illustrated in Fig. A3. As marching proceeds along the axial direction, the extent of the 
wake cut grows within that patch. The nomenclature for the starting and ending points of 
the wake cut are also indicated in Fig. A3. The number of points in the patch containing 
the wake cut is not allowed to change during the calculation. Thus, while exercising the 
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respace option in the region containing the wake, the user has to ensure that the number 
of points in the wake patch (usually there are two wake patches; one corresponding to the 
upper cut and one for the lower cut) is not altered. 

The shape of the trailing edge is provided by the user using the update option. 

For the wing-body-vertical case of Fig. A2, a 3-patch initial region, a 6-patch center 
region, and an 8-patch final region was used. The patch definition for Region 1 is as indi- 
cated in Fig. A5. Zero length patches are asl permissible. Since the analysis is marching 
in nature, a complete geometry data set is not required to begin and partially process a 
problem. Appropriate use of restart solutions allows continuation of the analysis as new 
or modified geometry becomes available. 

The format for a typical station is shown below. The group of cards is repeated for 
each station of a region. The last point of each patch (except for the last patch of a station) 
should have the same coordinates as the first point of the next patch. 


Card No. Format 

Field 

Name 

Description 

A1 F15.6,I5 

1 

XI 

The x value (longitudinal) 




of this station. 


2 

ISCl 

The number of patches for this 


section. 1 < ISCl < 15. 

The group of cards A2 through A3 are repeated ISCl times. 


A2 


315 


1 ITH Patch number <15. 


2 IPT Number of points in this 

patch. 2 < IPT < 30. 
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3 ND Mesh spacing parameters*. 

Typically the same for all 
stations of a region. 

The A3 card is repeated IPT times 

A3 2F15.6 1 YK Vertical location of point 

(positive upwards). Points start 
at top centerline (see Fig. A5). 

2 ZK Spanwise location of point. 

Cubic spline interpolation is performed on input patch data to derive the geometry. 
Linear interpolation is performed to define the geometry at a marching plane between 
input stations. 

Sample geometry data for the problem of Fig. A2 is presented in Table 1 and was 
developed using CDS 18 . 

Update File Directives 

The SIMP code is not intended to compute all cases without the user having to interact 
with the code. There may be cases which will require the user to incorporate specific 
changes to the code to obtain a solution. The changes that are frequently encountered are 
listed here. 

1) Shape of the trailing edge for wake calculation. In order to initiate the wake calcula- 
tion, the code has to know the starting and ending K values of the wake region in the 
marching plane (see Fig. A3). Depending on whether the trailing edge has a positive 
slope or a negative slope, the K values in Fig. A3 are properly computed by the code. 
The user has to prescribe the shape of the trailing edge. This update change is done 
in Program MAIN by prescribing Zte — /(?) = of + b where “a” defines the trailing 

* For segment AB: 0 equal space; 1 cluster near A; 2 cluster near B. 
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edge slope. The location where this change is made in MAIN is clearly marked in the 
code. If the trailing edge has no sweep, then KWKEDl = KWKSTl = K of leading 
edge (KLE). 

2) Respecification of the relaxation parameter OP in the grid generation routine. For 
severe geometry cases, to impose orthogonality and required specific grid spacing 
in the T) direction will require smaller values of OP (as low as 0.005). Generally, OP = 
0.02. The parameter OP provides underrelaxation for the constraints P and Q which 
impose the orthogonality and grid spacing near the body surface. When OP = 0, 
the grid solver does not impose any constraints. For geometry sections with drastic 
variation in slopes in the cross section, imposition of P and Q constraints can lead to 
instability in the grid solver. Lowering the value of OP should relieve this problem. 

Usually, the grid in a marching plane is divided into several subregions (the NRM 
parameter defines the number of subregions, see Fig. A3). The value of OP can be set 
to different values for different grid regions depending on the severity of the geometry 
slope change in that particular region. 

The update change for OP is made in Subroutine GRID in loop “DO 100”. 

3) Averaging of In regions of rapid flow expansion (flow around a sharp leading edge, 
clipped wing tip) the marching algorithm can run into instability problems. This will 
result in the density value at certain grid points becoming negative. Even though 
this might occur at only one or two points at a marching plane, the whole marching 
calculation comes to a halt. This problem, for most cases, can be alleviated by reset- 
ting the potential <f> values, at grid points encountering negative density, by averaging 
the <f> values surrounding those grid points. One such <f> averaging technique usually 
specified through the update option is illustrated in Fig. A7 for a sharp leading edge 
point. If the problem persists, even after implementing the <f> averaging procedure, 
one might consider a similar averaging procedure for density also. The user should 
carefully examine the grid to make sure that the negative density problem is not due 
to improper grid or improper geometry. 
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The update for <f> averaging is made in Subroutine INVSETA. The location is after 
the statement “705 CONTINUE”. 


4) Restart for nacelle on. In order to introduce a flow through nacelle in the calculation, 
the respace/restart option is exercised just at the nacelle face location as illustrated 
in Fig. A8. The calculation is first performed just prior to the nacelle face location 
(Region N). Then, in the overlapped region, the nacelle face is extended forward 
(fictitious flow through nacelle surface) and a respace is performed to set up grid and 
<f> values that correspond to the nacelle-on geometry. Such a respace calculation with 
a fictitious flow through nacelle is essential for a smooth transition into nacelle-on 
marching calculation. Through update, the user has to prescribe the K value of the 
starting and ending points of the nacelle to aid in the respace procedure. A sample 
update procedure for a nacelle-on calculation is given in Appendix B for the Zone 3 
calculation. 
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0000 - 
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0240= 

02503 

0260* 

0270= 

0280= 

02003 

0300* 

0310= 

0320= 

0330= 

0340= 

0350= 

0360= 

O370« 

0300= 

0300= 

0400* 

04103 

0420= 

0430= 

0440* 

0450* 

0460- 

0470* 

0480* 

0400- 


1.424664 
1.356820 
1.204987 
1.005808 
0.017631 
0.640761 
8.257428 
-0,130001 
-0.450283 
-0.64S5S0 
-0.004200 
-1.161518 
-1.204488 
-1.417318 
-1.456881 
-1.600636 
-1.781206 
-1.014585 
-1.060516 
2S7. 762030 

16 0 

60.010070 

50.010070 
40.824403 
40.187645 
40.0331 IQ 
46.350802 

2 6 8 

46.350802 
44.232117 
41.780586 
80.281430 
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4.056256 

7.001464 

18.603526 

13.506638 

16.814454 

16.814454 
18.186864 
10.800887 
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Fig. Al. Flow chart for the full potential code. 
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KWKST = KWKMIN; KWKED = KWKMAX 

K VALUES OF KWKEDI AND KWKSTI ARE COMPUTED FROM ZWAKEI 



Fig. A3. Cross-section patches and nomenclature. 
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NEARLY CONICAL SHAPE 
CALCULATION BASED ON CFLNO 



NONCONICAL SHAPE 

Af BASED ON DZMIN AND DZMAX 


Fig. A4. Marching step size selection. 
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Patch 3 


Fig. A5. REGIONl patching. 
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REGION 2 


Fig. A6. Cross section patches in overlap region. 
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<p[ K, 2) = [0(K - 1, 2) + 0(K + 1,2) + <P( K, 3)]/3. 


Fig. A7. Potential averaging at wing tip clip region. 
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WITH FICTITIOUS 
NACELLE EXTENSION 
(FLOW-THROUGH 
NACELLE LINE) 


Fig. A8. Nacelle-on calculation. 


APPENDIX B — DEMONSTRATION OF A TEST CASE 


A sample test case is presented here to familiarize the user with the operation of the 
code. 

Figure Bl shows the surface grid for the fighter configuration discussed in Appendix A. 
The configuration is an advanced fighter concept consisting of a blended wing/body with 
underslung nacelles and twin vertical tails. As mentioned in Appendix A, the solution 
for a complex configuration is accomplished by breaking the configuration up into regions 
or zones. These separate calculations are necessary because the configuration becomes 
increasingly complex in cross-sectional shape as marching is done in the axial direction. 
To accommodate the emergence of the wing, nacelle, tail, and wake region requires that the 
number of points in the circumferential direction be increased from zone to zone through 
use of a respace option. 

A brief description of the appropriate header information for each analysis zone is 
given below. The test case is for a combined yaw and angle of attack flight condition, 
M = 1.6, (3 = 4°, and a = 6°. 
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Zone 1 — Input Description 

This zone initially performs the conical calculation necessary to set up the starting 
plane solution needed to initiate the marching calculation. The parameter NCON specifies 
the number of conical iterations to be performed for the crossflow geometry at ZTAl. 

The header data for this region is given and a display of the flow field pressure contours 
and grid at the end of Zone 1 is shown. A sample of the printed output data and a discussion 
of the display postprocessor will be given at the end of Appendix B. 



Cross-section patch definition for Zone 1. 
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Zone 1 — Flow field pressure contours and gridding at x = 370. 


Zone 2 — Input Description 


This zone is required to increase the number of circumferential points to accommodate 
definition of the wing. The number of cross-sectional patches increases from three to six 
and there is a slight overlap region to adjust the Zone 1 flow field solution to the new grid. 
The header data is given along with a flow field display at the end of Zone 2. 



Cross-section patch definition for Zone 2. 
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Zone 2 — Flow field pressure contours and gridding at x = 450. 


Zone 3 — Input Description 


Zone 3 is required because the underslung nacelles on the configuration will be mod- 
eled as flow-through ducts. The primary changes required other than the obvious one 
of modifying the contour geometry is to overlap Zone 2 and Zone 3 both upstream and 
downstream of the inlet face. This permits both transition between six and nine patch 
definition required to add the two sides and the bottom of the rectangular-like nacelle, and 
allows associated grid index information and temporary code corrections to be supplied 
via the update file. 

Again, the header data is given as well as the necessary update statements and a flow 
field display at the face of the nacelle. 



Cross-section patch definition for Zone 3. 
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Zone 4 — Input Description 


Zone 4 is required to incorporate the vertical tail and the wing wake. This region 
corresponds to the region 3 solution 2 discussion of Appendix A. The number of patches 
is increased to 11 and the solution of Zone 3 is interpolated on to the new respaced grid 
through an overlap region between Zone 3 and Zone 4 . 

The total solution is now composed of four sequential regional solutions for Zones 1 
through 4 . The header data for Zone 4 is given with a flow field display approximately 
midway in the region. 



Cross-sectional patch definition for Zone 4 . 
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Zone 4 — Flow field 


pressure contours and gridding at x = 680. 


Output Data Description 


Sample printed output data is presented on Table Bl. Standard tabulated data is 
produced every NP marching steps as defined in the header data. More detailed physical 
plane data can be output at specified stations using the parameter NSPTI and the tabu- 
lated stations following the true/false header data input. Cartesian coordinates, velocities, 
and mesh indices are indicated in the following sketch. The axial velocity component, U, 
is positive out of the plane of the paper. 


A v V 
I X ' V 

k=2 



Crossplane data and tabulated surface pressure coefficient data at constant span sta- 
tions (i.e., z) are output via Tape 1. These files are displayed using the postprocessor 
described in the next section. 
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TABLE Bl CONCLUDED 


Postprocessor Description 


A Tektronix based graphics display capability is available for displaying results. The 
file containing the cross-sectional data is generated for ZTA1 by the supersonic full poten- 
tial analysis as Tape 1. 

Typical program prompt response is 


X-STA7I0N • 54 .55 


IPL0T...I5 

• 2 FOR MACH NO. CONTOURS 

• 3 FOR PRESSURE CONTOURS 

• 10 FOR CCHPUTATIONAL GRID 

• 12 FOR HACh NO. PROFILES 

• 13 FOR PRESSURE PROFILES 

e — + — -e — o — * — e — 0 

’ le 

xniN.xnftx.vniN.YMPx. . .4Fio.<* 

X is THE PLOTCS X COORDINATE 

Y IS THE PLOTCS V COORDINATE 

V IS THUS MACH NO. FOR 1PL0T-12 

Y IS THIS PRESSURE FOR IPLCT-12 

V IS OTHERL’ISE THE PHYSICAL Y COORDINATE 

X IS ALUAYS THE PHVSICAL X COORDINATE 
xnN. 0. Xf-lAX* .22TS70S3E+03 

VMIN* - .22351750E a 03VCAX* .2293S620E-t03 

0 4 . 0 4 0 4 0 4 0 

7 0.0 25. -10. 10. 


Typical output is shown below. The results are truncated at Z = 25 and Y = ±10 as 
specified in the input above. 


n*CH no. • l.U 

«IPMA • (.11 


CO^lTbtKna; c*i; 
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HACH NO. 
ALPHA • 


X-5TAT10N • 54. 9S 

IPL0T...I5 

• 2 FOR MACH NO. CONTOLRS 

- 3 FOR FRESSuRE CONTOLRS 

- 1C FOR COMPUTATIONAL GRID 

- 12 FCR HACH NO. PROFILES 

• 13 FCR PRESSURE PROFILES 

0 * 0 + 0 ♦ 0 ♦ 0 

’ 3 

xniN,xriAx,Ynih,Yr;AX. . ,4Fie.4 
X IS THE PLOTCS X COORDINATE 

V IS The PLOTCS Y COORDINATE 

V IS Thus HACH NO. FOR IPL0“-12 

Y IS THUS PRESSURE FOR IPL0T-13 

Y IS OVERUSE THE PHYSICAL Y COORDINATE 

X IS ALUAYS THE PHYSICAL X COORDINATE 
xniN* 0. XHAX » .22787093E+03 

VPIN- - .22351 750E* 03 YflAX* .22933620Et03 

0 + 0 ♦ 0 ♦ 0 ♦ 0 

•9 0 .C 25. - 10 . ie. 

FniN- -.74e3S2iE-ee fhax- . 2501 E 63 E +00 
CniN,CnAX # DCL...3 c 10.4 
CHIN* LCUEST CONTOUR LEjEL 
CP1AX* HIGHEST CONTOUR LEUEL 
DCL ■ INCRETENT BETUEEN CONTOUR LEVELS 
0 0 + 0 4 0 «► 0 

* -.£4 .25 .es 

SYMBOLS ARE TO BE PLOTTED ON EJERY (12) CONTOUR LINE. 

* 02 

SYMBOLS ARE TO BE PLOTTED EUERY (12) POINT (5 ) . 

* 02 


- 1.28 
5.00 


PRESSURE CONTOURS 


CCntouR LEVELS 
LEVEL INC • 0.050 



10 


15 


£0 


25 
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X-STATION • 54. S5 


1PIOT ... 15 

• 2 FOR RACH NO. CONTOURS 

• 3 FOR PRESSURE CONTOURS 

• 10 FOR COMPUTATIONAL GRID 

• 12 FOR RACH NO. PROFILES 

• 13 FOR PRESSURE PROFILES 

0 — — 0 1 0 ♦ 0 — 0 

? 13 

FRIN* -,7403621E*O0 FRAX- .25018S3E+00 
XniN.XRAX, VfllN.YRAX. . .4F10.4 
X IS THE PLOTtS X COORDINATE 

Y IS THE PLOTCS V COORDINATE 

Y IS THUS HACH NO. FOR IPLOT-12 

Y IS THUS PRESSURE FOR IPLOT-13 

Y IS OTHERUISE THE PHYSICAL Y COORDINATE 

X IS ALUAVS THE PHYSICAL X COORDINATE 
XJ1IN* 0. XRAX* .22787093E+03 

YRIN* -.223S1750E+O3YRAX* . 22938G20E+03 

0 ♦ 0 ♦ 0 ♦ 0 — — 0 

7 0.0 25. 
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Conservative Full Potential, Implicit Marching Scheme 

for Supersonic Flows 


Vijaya Shankar* 
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An aerodynamic prediction technique based on the full potential equation in conservation form is developed 
for the treatment of supersonic flows. This technique bridges the gap between simplistic linear theory methods 
and complex Euler solvers. A local density linearization concept and a second-order-accurate retarded density 
scheme, both producing the correct artificial viscosity, are introduced in developing an implicit marching scheme 
for solving the scalar potential </». Results for conical flows over delta wings, cones, and wing-body com- 
binations, and for nonconical flows over bodies of revolution at angles of attack are compared with Euler and 
nonconservative full potential calculations and experimental data. The present formulation requires an order of 
magnitude less computer time and significantly less computer memory over Euler methods. 


I. Introduction 

A ERODYNAMIC prediction techniques that can handle 
significant geometric complexity for use in supersonic or 
hypersonic configuration design are based on either hyper- 
sonic impact methods 1 or linear theory analysis, 2 both of 
which require minimum response time and cost. However, 
shortcomings are present in both the impact and linearized 
methods. Aside from these simplified techniques, limited 
capabilities also exist for calculating supersonic flowfields 
using very complex Euler codes, 3 ' 6 using either shock cap- 
turing 3 or shock fitting 4 ' 6 methods. The use of these codes as 
viable aerodynamic prediction techniques for configuration 
design is, however, not practical due to their slow response 
time (requirement of large computer memory) and excessive 
computer cost per run due to strict stability requirements. 
Thus, we have on one end of the spectrum, very simplified 
codes that require minimum computer time to provide less 
accurate results and, on the other end, very complex Euler 
codes that require excessive computer time to provide quality 
results. 

In an attempt to bridge this gap between simplistic linear 
theory methods and complex Euler solvers, several 
methodologies such as the second-order potential analysis, 7 
hypersonic small disturbance theory, 8 and, more recently, 
nonconservative full potential methods 9,10 have been con- 
sidered by various investigators. The second-order theory, 7 in 
spite of the significant improvements reported, suffers from 
the lack of nonlinearity in resolving proper cross flow shocks 
and sonic lines. Also, the singularities inherent in the for- 
mulation create difficulties in the numerical treatment of 
subsonic leading edges. The finite difference analysis of the 
hypersonic small disturbance theory 8 indicates that the 
solution procedure is as complex as that for the Euler 
equation and not particularly responsive to preliminary design 
level of effort. 

Recently, Grossman 9 and Grossman and Siclari 10 have 
computed supersonic flowfields over conical and nonconical 
cambered and twisted delta wings with remarkable success 
using the nonconservative full potential equation and a 
transonic relaxation method. However, their approach is 
made complicated by the use of global conformal mappings 
which apply only to certain classes of configurations. Also, 
the nonconservative form of the full potential equation is in 
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terms of second derivatives of the potential <t>, which, when a 
transformation is applied, generates a large number of first 
and second derivative transformation terms. 

The full potential method proposed in this paper is 
significantly different from that of Refs. 9 and 10. First of all, 
the method is based on the conservative form of the full 
potential equation, since for a shock capturing procedure to 
conserve mass across the shock wave, 11 it is essential that the 
equation be cast in conservation form. 12 Second, the method 
can accommodate a numerical or analytical mapping 
procedure that is either orthogonal or nonorthogonal without 
complicating the form of the equation, in contrast to Refs. 9 
and 10. Third, the method is based on an approximate fac- 
torization implicit algorithm that can yield convergence much 
faster than the conventional successive line over-relaxation 
method. 10 Finally, the method is not an adaptation of a 
transonic code using type dependent operators, but a scheme 
specifically developed and tailored for supersonic marching 
problems using a density linearization concept and has no step 
size restrictions. 

To validate the present methodology, results are shown for 
a variety of conical and nonconical geometries and are 
compared with Euler solutions and full potential results of 
Refs, 9 and 10. Results indicate that the method works just as 
fast and efficient for nonconical flows as in the case of conical 
geometry treatment. Results also indicate that the method is 
very useful in computing very high-speed flows (M v ^ 2-6) for 
the moderate flow deflection angles (a -4-10 deg) where the 
neglect of entropy generation does not seriously distort the 
main features of the flowfield. 

The present method can also handle more complicated 
geometries (realistic wing-body combinations) than the ones 
reported in the paper, but requires a suitable grid generation 
routine, especially near wing-body junction regions. In a 
subsequent paper, 13 results for nonconical wing-body flows 
will be presented along with a formal method of charac- 
teristics treatment for cross flow signal propagation. 

II. Formulation 

The conservative form of the full potential equation in 

Cartesian coordinates x , y, z can be written as 

d(pu) t d(pt>) | d(pw) _ G (1) 

dx dy dz 

where p is the density and w, v , w are the velocity components. 

They are calculated as the gradient of the potential <£, 

u = <t> K .\ v = <t> v ; 


w = <£. 


( 2 ) 
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The density p is computed from the isentropic formula 

/ — Ml>(u 2 + v 2 + w 2 - 1)\ (3) 

If the density is normalized with respect to the freestream 
value, then the speed of sound a is given by 

(4) 

where A/* is the freestream Mach number. 

The objective of the paper is to solve for the scalar potential 
<t> from Eq. (1) subject to the surface tangency condition 
<t>„ =0 (n is normal to the body surface). Examining Eq. (1), it 
is very clear that <f> appears in a nonlinear form due to the 
presence of the density term inside the derivative. The ap- 
proach to be described here is a method that treats the density 
term in such a way that it produces the correct artificial 
viscosity needed for shock capturing and that enables one to 
solve for <t> with relative ease. 

In order to apply the surface tangency condition at the 
actual body location, a body-fitted coordinate transformation 
is essential. Introducing a body-fitted coordinate trans- 
formation, f=r(jr,j\;:), rj = ri(.x,y,z)* and £ = £(*,>%;:), Eq. 
(1) transforms to 

(«j),4D,47 ),-• 

where U, V % and W are the contravariant velocity com- 
ponents. Introducing the following notation for convenience: 



s 

II 

w=u s 

II 

il 

* 

II 

li 

<3 

II 

* 

z=x 3 


the contravariant velocities and density are given by 

3 

u,= E <*,,**, 

j-i 

y dX, 3JT 
a, ‘~h Sx„ dx„ 


;= 1.2.3 

/= 1,2,3, j= 1,2,3 


The Jacobian of the transformation J is represented by 


which indicates that for stability, the form of artificial 
viscosity be 

U *m + + 1 

+ A{ W[ U<t> {H + )) (8) 


assuming that U, V, W are positive. What this implies is that if 
the flowfield is hyperbolic (q>a), then solution can be ob- 
tained by marching along the hyperbolic flow direction s. 
Once the total velocity q becomes less than a , then marching 
along s is not possible. This is reflected in the fact that the 
effective artificial viscosity given by Eq. (8) is now negative. 

Now we will proceed with the numerical procedure for 
solving Eq. (5), and show the resemblance of the resulting 
artificial viscosity to that of Eq. (8). 


A. Treatment of inEq.(5) 

*r 

Consider the direction fto be the marching direction. The 
condition to be satisfied for this to be true will become evident 
at the end of this analysis. Both the density p and the con- 
travariant velocity U are functions of the potential 4) and the 
transformation metrics, as represented in Eq. (6). In order to 
finite difference this f derivative quantity in terms of 4> only 
will require some linearization treatment of the density. This 
will be termed the “local density linearization” procedure. In 
the transonic formulation described by Holst, 15 the density is 
upwind biased and computed at the old level, while retaining 
central differencing for the (U/J) term at the current level. 
Such an upwind density bias is shown to produce the right 
artificial viscosity in Ref. 14. Referring to Fig. 1, for a pure 
supersonic marching problem (say we want to march from the 
ith plane to the / H- 1th plane), a transonic relaxation 
procedure 15 in the marching direction f is not appropriate 
because the solution <t> at the / + Ith plane is not influenced by 
the / H- 2th plane. Hence, the following marching procedure is 
developed. 

Given the <j> information at all the previous planes /, / — 1 , 

i-2 the problem is to compute <j> at the / -h 1th plane. 

Now, expand the unknown p=p(4>) in terms of a neighboring 
known state denoted here by a subscript 0 (ith plane in- 
formation would represent the neighboring known state for 
the /+ 1th plane). 


P = Po + 



£<t> + ... 


(9) 


3(r>i?.o 

d(x.y,z) 


r, r, ' 

Vx Vy Vz 

tr Zz . 


(7) 


Equation (5) is in terms of a general coordinate system 
(T* 1 ?^) an d can accommodate any kind of mapping 
procedure, either analytical (conformal mapping) or 
numerical type. Any numerical marching procedure applied 
to Eq. (5) to simulate a supersonic flow should have a 
truncation error whose leading terms represent a correct 
artificial viscosity. This is essential to ensure marching 
numerical stability and to exclude the formation of expansion 
shocks which are unphysical and correspond to a decrease in 
entropy. The nature of the required artificial viscosity can be 
studied by an analysis 14 of the canonical form of Eq. (5), 


KNOWN DATA 



Fig. 1 Implicit computational molecule. 
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where A0 = 0 - 0 O and dp/d<t> can be shown to be a differential 
operator 16 


( — ) A<t>=- P - f 

r a a 

+ W 0 3z\ 

\d<t>/ 0 a 2 0 1 

l 0 ar dr. 

0 as J 


o) (10) 


Substituting Eqs. (9) and (10) into the first term in Eq. (5), we 
get 


djpU/J) d p 0 

ar ± afvr < ’ aj, 



( 11 ) 


Substituting for U in terms of 0 from Eq. 6 and rearranging 
Eq. (1 1) in terms of the potential difference A0, we get 


f Po 

7 u ! ' 

\ a / 

UV\ 

It! 





ar 

d / UW\ d l') 

X-A0+(a„-— ) -A* + t/ 0 ]) 


(12) 


where the speed of sound a 0 , the density p 0 , and the con- 
travariant velocities U 0t V 0 , W 0% represent information at the 
neighboring known plane. The f derivative term of Eq. (12) 
will now be one-sided differenced. Assuming t/ is positive. 


=f r i i (i3) 

An upwind differencing of the form Eq. (13) applied to Eq. 
( 1 2) can be shown to produce a truncation term whose leading 
term is 



which will always represent a positive artificial viscosity as 
long as 

U 2 

— >a 2 (15) 

a u 

The preceding relation sets the condition for f to be the 
marching direction, for if U 2 /a n is less than the square of the 
local speed of sound, then the artificial viscosity becomes 
negative and a marching instability will occur. This also 
implies that the projection of the total velocity vector q in the 
direction normal to the const plane (t;,£ plane) is super- 
sonic. For example, in a spherical system, for the radial 
direction r to be the marching direction, the radial velocity q r 
must be supersonic. The similarity between the artificial 
viscosity term given by Eq. (14) and the first term appearing in 
Eq. (8) can be noted. When backward differenced, the terms 
in Eq. (12) will lead to a diagonally dominant tridiagonal set 
of equations for the unknown A0 when coupled with the other 
two terms in Eq. (5). The mixed derivative terms like <f>^ and 
<t>x appearing in Eq. (12) will be upwind biased, depending on 
the sign of the coefficient multiplying them to preserve 
diagonal dominance and to provide the right artificial 
viscosity. 


„ d(pV/J) 

B. Treatment of in Eq. (5) 

dv 

This term will be written at the /+ 1th plane to make the 
resulting scheme fully implicit. 


The density term p in Eq. (16) cannot be represented at the 
/H- 1th plane since that would result in a very complicated 
nonlinear form for 0. Hence, a density approximation is 
introduced by setting p=p* where p m -p 0 for conical flow 
treatment. In the case of nonconical flows, while advancing 
from /' to the /+lth plane, several iterations are performed 
within each cross flow plane (77, £) to refine the density p * to 
properly account for the axial geometry variation. This is 
done by initially setting p* to p 0 and then subsequently 
refining it by setting p* to the previous iterate value of p at the 
current plane. In many cases where the axial variation of the 
geometry is gradual (especially for smaller step size 
calculations) it was found that setting p* =p 0 even for non- 
conical flows produced very good results without having to 
refine the density subsequently. 

Writing Eq. (16) in terms of the potential difference A0 


d{p V l J) / p*a 2 , A 0\ / p*a 22 d \ 

dr, ~\ J A f/, V J dr, V, 


/p*a v d \ (P*Qt> d \ /P*a n d \ 

(Ti“)*(TE‘').*(Tii‘'). < 17 ’ 


A simple central differencing for the various terms in Eq. 
(17) will not be sufficient as that would not provide the 
desired artificial viscosity given by Eq. (8), required for shock 
capturing. To simulate an artificial viscosity of the form given 
by Eq. (8), the density will be upwind biased based on the 
previous work reported in Refs. 14-16. The density p* will be 
replaced by a modified density p* given by 


+ Mrj+'-'.kl </ + *> <pV-’»u + U-*) <P*)y-/ + *M ) U8) 

where m - 0 when (K 0 ) y+ k .>0 and w= + l when 
(Vo) j+'' : ,k <0. When 6 is set to zero, first-order accurate 
density biasing is achieved while 9 = 2 gives second-order 
accuracy. The artificial viscosity coefficient y j+ ^ k is com- 
puted as follows: 

+ U ( ,9 > 

where 5 = 0 for V j¥ ,. s tk >0 and s= 1 for ,,^ k <0. 

Treatment of density as represented by Eqs. (18) and (19) 
would always produce a positive artificial viscosity as long as. 
the local total velocity q 0 is supersonic. If that becomes 
subsonic, then the marching procedure would fail and the 
problem have to be treated as a transonic problem. 

The treatment of the {d/d%)[pW/J\ term in Eq. (5) is very 
similar to the just described (d/dr,) [p V!J] term, except that 
the density biasing will be in the £ direction and will be based 
on the sign of W. 

C. Implicit Factorization Algorithm 
Combining the various terms in Eqs. (12) and (17), and the 
terms arising from (9/d£) (p W/J\ will result in a fully implicit 
representation of Eq. (5) which cannot be solved without 
introducing an approximate factorization procedure. After 
some rearrangement of the terms, the factored implicit 
scheme becomes 

r / 3 Fjhi j i 

L 0 d$ \ J Af / 0 J as J 

x r /+ dL 1 + L 1 (?£») 

L 0A$* dr, (3 dr, V JA$ / 


d(p V/J) 
dr, 


d_ 

dr, 


( a :i < i>t+ a n < t , T, + a 23 < l > c)} 


06) 


/ d p*a 
(3 dr, J 


22 _ 


d_ 

dr, 


} 


A d) = R 


(20) 
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This equation has the form 

L { L 1 (A0) =R (21) 

and it is implemented as follows: 

L £ (A0)*=/? L,(A0) = (A0)* 0 = 0„ + A0 (22) * 

The various quantities appearing in Eq. (20) are given by 



and the right-hand side term R consists of various known 
quantities. 

The algorithm Eq. (22) requires only scalar tridiagonal 
inversions. Also, the scheme does not pose any restrictions on 
the direction of sweep that are present in the successive line 
over-relaxation method. 9 ,0 


D. Freestream Truncation Errors 
To subtract out any numerical truncation error due to 
incomplete metric cancellation, 16 it is essential to add the 
terms (especially for a highly stretched nonorthogonal grid) 


a /p.i/.- 

u 3 ( p - y -' 


ar v j > 

' drj \ J > 

> + d£ V J 


) 


(24) 


to the right-hand side of Eq. (20). 


E. Boundary Conditions 

In order to solve for A0 from Eq. (20), boundary con- 
ditions will have to be prescribed at all four boundaries as 
shown in Fig. 2 at the current /+ 1th plane. While performing 
the Z. £ operator in Eq. (21), boundary conditions in terms of 
A0" will be required along the k = 2 and *=AMAX-l 
boundaries. For a pure angle-of-attack problem, k- 2 and 
£ = /CMAX-l can be considered as planes of symmetry 
across which all flow variables reflect. The quantity A 0*, even 
though it has no physical significance, can be safely set 

( A0 ) * + IJ'KMW = (^0 ) U Ij'.XMAX - 2 

(A0)f +/i ,./=(A 4>)U U j (25) 

Jhe L n operator would require boundary conditions along 
j -2 and y=/MAX in terms of A0. Since y = 2 is the body, the 
surface tangency condition 

K=tf 2/ 0 f + ff; 2 0,+ffjj0|=0 (26) 



j 2 IS BODY 

, * l.k = 1. k - KMAX ARE DUMMY ARRAYS 
I = JMAX IS FREE STREAM 

Fig, 2 Physical and computational plane. 


will be set at all points (1+1,2,*). Along j = J MAX, 
freestream A0 will be imposed. 

F. Grid System 

As shown in Fig. 2, the physical space ( x,y,z ) is trans- 
formed into a body-fitted (f,ij,£) computational space. The 
transformation is performed numerically by using the elliptic 
grid generation techniques originally developed by Thompson 
et al., 17 and later modified by Steger and Sorenson, 18 and 
Middlecoff and Thomas. 19 The present full potential method 
does not require an orthogonal grid; however, the error in- 
troduced by the approximate factorization, Eq. (20), can be 
minimized if the grid is orthogonal in the cross flow plane 
(* 7 ,£) . For conical flow calculations, the grid is generated only 
once, and, as the marching procedure continues, the grid is 
allowed to grow conically. For a general nonconical body, it 
would be necessary to construct the grid in every marching 
plane. 


III. Results 

Results are presented for both conical and nonconical 
supersonic flows. Comparisons are made with Euler 3 - 3 and 
full potential 9,10 results and experimental data. All the 
calculations were performed using a CDC 7600 machine. 

A. Conical Flows 

Besides validating the methodology, computation of 
conical flows is of interest for generating the initial data plane 
for nonconical calculations. For a conical geometry (radially 
invariant), the initial data plane with freestream conditions is 
chosen at some location f=f 0 (usually set at $*=1). The 
solution is then marched along fusing Eq. (20) and boundary 
conditions. The conical flow calculation is assumed to have 
converged when the change in the root mean square density is 
less than 10 ~ 4 . 

Supersonic Leading- Edge Delta Wing 

Figure 3 shows the compression surface pressures for a 
supersonic leading-edge delta wing at =6, angle of attack 
- 8 deg, and leading-edge sweep A = 70 deg. The present 
full potential solution compares well with the Euler solution 19 
and experimental data. Also shown are the results from the 
first- and second-order linear theory. 7 Using a 30x40 grid in 
the (tj,£) plane, the present approach required 40 iterations to 
achieve convergence, and 12-15 s of computer time to produce 
the results shown in Fig. 3. 

It is interesting to note that in spite of the limitations of the 
full potential theory, even at a very high Mach number of 6, 
the comparison is in reasonable agreement with the Euler 



Fig. 3 Predicted compression surface pressure for a 70-deg sweep 
delta wing at = 6, a = - 8 deg. 
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solution and is significantly better than the second-order 
theory. The discrepancy between the full potential and Euler 
results is mainly due to neglect of entropy generation in the 
present approach. 

Circular Cone and Ellipse 

Figure 4 shows the surface pressure distribution for a 
circular cone, half angle 7.5 deg, A/* =3 at 15 deg angle of 
attack. At this angle of attack the cross flow Mach number 
becomes supersonic as the flow turns around the cone from 
the windward symmetry to the leeward symmetry. This cross 
flow supersonic region is terminated by the formation of an 
embedded shock on the cone surface. This is evident from the 
results of Fig. 4. Grid clustering near the cross flow shock was 
used both in the Euler calculation of Kutler, 3 and in the 
present method, to finely resolve the pressure jump. The 
present calculation required 25 s of computer time, using a 
30 x 60 grid in the (17,$ ) plane. 

The liftoff of the vortical singularity on the leeward 
symmetry plane associated with the formation of the em- 
bedded shock is shown in Fig. 5. The location where the 
contravariant velocity V goes through zero on the leeward 
symmetry plane (1K=0) denotes the location of the vortical 
singularity. The behavior of the cross flow streamlines 
converging to the vortical singularity is also shown in Fig. 5. 

Figure 6 shows the full potential and Euler cross flow Mach 
number contours for the circular cone case. The presence of 
the embedded shock wave in both the results is very clear. The 
location of the vortical singularity liftoff is also shown in the 
figure. The Euler result is very oscillatory near the vortical 
singularity location while the present method predicts a 
smoother flowfield in the vicinity of the vortical singularity. 

The surface pressure distribution on an elliptic cone 
d c - 18.39 deg, 6 C = 3.17 deg at = 1.97 and a = 10 deg is 
shown in Fig. 7. The results of the present study are compared 
with Euler calculations of Siclari, 5 full potential results of 
Grossman, 9 and the linearized thin wing solution of Jones 
and Cohen. 20 The agreement between the various nonlinear 
methods is very good, including the position and strength of 
the embedded shock wave. 

Wing- Body Combination 

Figure 8 shows the numerically generated grid distribution 
in the cross flow plane (77,$) of a conical wing-body com- 
bination. The design of this conically cambered delta wing to 
achieve shockless recompression is reported in Ref. 21. Figure 
9 shows the pressure distribution around this wing-body 
combination at A/„=2, and a = 7.81 deg. The leading-edge 
sweep A is 57 deg. The comparison of the results from the 
present method with the experimental data 21 is excellent. The 
calculation used a 15 x49 grid in the (17,$) plane and required 
less than a minute of computer time. 




Fig. 5 Vortical singularity liftoff for a circular cone at \f 0B =3, 
0 = 15 deg, 6 C =7.5 deg. 



PRESENT EULER 

FULL RESULT 

POTENTIAL (REF. 3) 

RESULT 


Fig. 6 Comparison of cross 
flow Mach number contours 
for cone at angle of attack; 
Af* = 3, a =15 deg, B c =7.5 
deg (VS = vortical singularity). 



Fig. 4 Surface-pressure distribution for cone at angle of attack; Fig. 7 Surface-pressure distribution on an elliptic cone; A/ w = 1.97, 
M m = 3, a = 15 deg, 0 f = 7.5 deg. 6 C = 18.39 deg, h c = 3. 17 deg, a = 10 deg. 
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Fig. 8 Computational grid around a wing-body combination. 



Fig. 9 Surface-pressure distribution for a conically cambered wing- 
body combination; M „ =2, a = 7.81 deg, A = 57 deg. 


B. Nonconical Flows 

Results are also presented for nonconical bodies of 
revolution and compared with experimental data. The initial 
data plane for the nonconical marching calculation is first 
obtained by performing a conical calculation over an assumed 
very small conical nose. This conical calculation usually takes 
20-30 iterations on a typical 30x30 grid in the (jj,£) plane. 
The nonconical calculations did not exhibit any increase in 
computational time over the conical procedure. As mentioned 
earlier, for the applications considered here where the cross 
flow station does not vary substantially from the previous 
one, it was found that there was no need to iterate the solution 
at each cross flow plane (r j,£) to converge the density, and 
plottable accuracy was achieved by simply marching right 
along the body. However, if the body changes shape ap- 
preciably, the current implicit procedure might take 3-5 
iterations per cross flow plane to refine the solution. 

Reference 22 contains experimental data for several bodies 
of revolution at various Mach numbers and angles of attack. 
The shape chosen for comparison here is a circular arc- 
cylinder body. After the initial data plane was computed using 
a conical nose assumption, the current method typically used 
60 marching steps to reach the end of the body but the 
calculations are not subject to any step size restriction. A 
typical calculation required 40-45 s of computer time. Figure 
10 shows the circumferential surface pressure distribution at 
two different axial stations (x/f= 0.225 and 0.425) for 4 and 8 
deg angles of attack at Af*. = 2.3. The results from the present 
method are compared with experimental data, 22 showing very 
good agreement for the windward region with some 



Fig. 10 Circumferential pressure distribution for a circular-arc- 
cylinder body at Af* = 2.3. 



Fig. 1 1 Surface-pressure distribution for a circular-arc-cylinder body 
at o « 8 deg, = 2.3. 


discrepancy on the leeward side, possibly due to boundary- 
layer buildup. 

Figure 11 shows the surface pressure distribution in the 
axial direction along the windward and leeward plane of 
symmetry, at = 2.3 and a - 8 deg. Again, results from the 
present method compare very well with the experimental 
data. 22 

IV. Conclusions 

An aerodynamic prediction technique based on the full 
potential equation in conservation form is developed for the 
treatment of supersonic flowfields. A local density 
linearization concept and a second-order accurate density 
biasing scheme are introduced in developing an implicit 
marching procedure. The method produces results that 
compare well with Euler- solvers, and requires an order of 
magnitude less computer time and significantly less computer 
memory over existing Euler codes for the cases presented in 
the paper. In a subsequent paper, 13 results for more com- 
plicated nonconical wing-body flows are presented, along 
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with a formal theory for the characteristic signal propagation 
in the cross flow plane. 
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A nonlinear aerodynamic prediction technique based on the full-potential equation in conservation form has 
been developed for the treatment of supersonic flows. The method uses the theory of characteristic signal 
propagation to accurately simulate the flow structure, which includes shock waves and mixed elliptic-hyperbolic 
crossflow. An Implicit approximate factorization scheme is employed to solve the finite-differenced equation. 
The necessary body-fitted grid system in every marching plane is generated numerically, using an elliptic grid 
solver. Results are shown for conical and nonconical wing-body combinations and compared with experimental 
data and Euler calculations. The method demonstrates an enormous savings in execution time and memory 
requirements over Euler methods. 


I. Introduction 

N ONLINEAR aerodynamic prediction techniques based 
on the Euler equations 1 * 3 and the full-potential 
equation 4 * 7 are steadily maturing into complex aerodynamic 
tools and becoming an attractive alternate approach to using 
the linearized panel methods. 8 Panel methods can handle very 
complicated geometries requiring minimal computer time to 
provide less accurate results, while the Euler solvers need 
expensive computer runs even for simple wing-body con- 
figurations. The full-potential methods 6 - 7 are a substitute for 
the Euler methods 1 - 3 to avoid the requirement of excessive 
computer time and memory allocation. While using a full- 
potential method for supersonic flows, one should be aware 
of the isentropic limitations of the theory. As a general rule, 
the full-potential theory is expected to perform well when the 
product of the Mach number and the characteristic flow 
deflection angle is less than 1 (Af5s 1). 

The full-potential method of Refs. 4-6 is based on the 
nonconservative form of the equation, while Ref. 7 and the 
present paper deal with the conservative form, to conserve 
mass across the shock. 9 - 10 In order to properly treat the 
supersonic flow structure, which includes shock waves and 
mixed elliptic-hyperbolic crossflow, the present method uses 
the theory of characteristic signal propagation based on the 
eigenvalue system of the full-potential equation. An ap- 
proximate factorization implicit scheme, which includes a 
density biasing procedure in the crossflow plane, is in- 
corporated to accelerate the computational efficiency. The 
density biasing procedure is activated by the eigenvalue 
system and properly takes into account the direction of the 
crossflow. The implicit approximate factorization scheme 
does not pose any restrictions on the direction of sweep that 
are present in the successive line overrelaxation method 
(SLOR). 4 * 6 

The full-potential as well as Euler methods require the 
application of boundary conditions at the actual body surface 
location. This, in general, necessitates the use of a body-fitted 
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Thermophysics, Fluids, Plasma and Heat Transfer Conference, St. 
Louis, Mo.. June 7-11, 1982; submitted July 7, 1982; revision received 
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Astronautics, Inc., 1982. All rights reserved. 
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coordinate system. In the present method, the equation is cast 
in a more general arbitrary coordinate system and the ap- 
propriate body-fitted grid is generated numerically, em- 
ploying an elliptic grid solver. 11 

The paper presents various results for conical and non- 
conical wing-body configurations and comparison is made 
with experimental data and Euler solution. The effect of the 
density biasing based on the characteristic signal propagation 
is demonstrated in terms of a sharper pressure profile across 
the shock wave. References 5 and 6 present excellent results at 
low supersonic Mach numbers, while Ref. 7 and the present 
paper demonstrate the capability of the conservative full- 
potential approach in handling even very high Mach number 
flows (Af* -4-6, a -0-8 deg). All of the calculations reported 
in this paper were performed using the CDC 7600 computer 
and clearly demonstrated an order-of-magnitude or more 
reduction in computer time over Euler methods. A typical 
nonconical wing-body calculation takes less than 2 min of 
execution time to produce results comparable with ex- 
perimental data. 


II. Formulation 

The conservative full-potential equation cast in an arbitrary 
coordinate system defined by t={{x,y,z), ri = T){x,y,z), and 
£ = £ (x,y,z) takes the form 

(*j ),4J)-W).- 

where U, V, and W are the contravariant velocity com- 
ponents. Introducing the following notation for convenience 


u=u„ 

ii 

5 

II 

x=x { , 

H 

II 

M 

II 

w* 

r=*/. 

v=x 2 . 



the contravariant velocities and density are given by 

3 
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y, dXi dXj^ 

hi dx i dx k 


i= 1,2,3 . 

j= 1,2,3 


p=[i-( 1 y-)MU U4>t + y<t>, + W'*, -■m]"'’"'’ ( 2 ) 


The Jacobian of the transformation J is represented by 


3 (x.y.z) 


■fr r, f t ' 

lx r i , i, 

S X Zy Z z . 


(3) 


Equation (1) is in terms of a general coordinate system 
and can accommodate any kind of mapping procedure, either 
analytical (conformal mapping) or numerical. Use of Eq. (1) 
to simulate the supersonic flow by marching in the f direction 
first requires the establishment that the equation is indeed 
hyperbolic with respect to the marching direction. The nature 
of Eq. (1) can be analyzed by studying its eigenvalue system. 
Combining the irrotationality condition in the (f,r/) and (££) 
plane and Eq. (1), one can write the following matrix 
equation: 

Aq c +Bq, + C<? ( = 0 (4) 

where 


A = 


B= 


C= 


), f 

0 

0 

-/ 

0 

( UJ){p\v ) #r 
0 


-1 



U/J) (pU)^ 
1 
0 

U/J)(pV )*, 
0 
0 

u/J) 

0 

0 


U/J)(pU)^ m 
0 
1 

u/J){pV) H ' 
0 
0 

U/J)(pW) H ’ 
0 
0 


The subscripts in Eq. (4) denote differentiation with respect to 
that variable. 

In order for Eq. (4) or Eq. (1) to be hyperbolic in the f 
direction, the following two conditions must be satisfied: 

1) A ~ l mu^t exist. 

2) All real linear combinations of and A~ l C must 
have real eigenvalues (characteristics). This implies A~ 1 
(aS+/3C) must have real eigenvalues for all combinations of 
a and & satisfying a 2 +& 2 - 1 . 

When the two conditions are applied to Eq. (4), the 
following criterion is obtained for f to be the marching 
direction. 

(pt/), f =p(<7„-^)<0 (5) 

where the transformation metric a n is defined in Eq. (2) and a 
is the local speed of sound. Equation (5) is the most general 
form. For example, in a spherical system (r,0,<£), for the 
radial direction r to be the marching direction, according to 
Eq. (5), the radial velocity q r must be supersonic. In a Car- 
tesian system (x,y,z), for x to be the marching direction, the 
velocity u must be supersonic. For convenience, the derivation 
of Eq. (5) for a Cartesian system is described in Appendix A, 
and the derivation for an arbitrary coordinate system (f,i 7 ,£) 
has been derived in a similar manner. 


Thus far, the condition for f to be a marching direction has 
been identified from the characteristic theory. This means the 
( 17 , £) plane will be treated as a marching plane, which will be 
defined from here on in this paper as the crossflow plane (the 
real crossflow is the projection of the velocity vector on a unit 
sphere with center at the origin). Even though the flow is 
supersonic in the marching direction (i.e., hyperbolic type), 
the behavior of the flow structure in the crossflow plane ( 77 , £) 
can be a mixed elliptic-hyperbolic type. Depending on the 
nature of the flow at a crossflow plane grid point (whether 
elliptic, parabolic, or hyperbolic), the 77 and $ derivative terms 
in Eq. (1) will be appropriately modeled. Again, the theory of 
characteristics will dictate how the signals are propagated in 
the crossflow plane. 

A. Crossflow Signal Propagation 

The nature of the flow in the ( 77 , £) plane can be analyzed by 
separately studying the eigenvalues of A ~ l B and A “'C. The 
eigenvalue character of A~ l B will determine the 77 -derivative 
treatment and similarly A~*C for the £ derivative. For 
illustration, only the study of A~ l B is shown here, and^l^C 
follows the same procedure. 

The eigenvalues X of A~*B are obtained by setting the 
determinant \A~ ! B-\I\ =0. Since A‘ / is assumed to exist 
[condition 1 preceding Eq. (5)], the following is true. 



CASE 1. ELLIPTIC 

CROSSFLOW 




CASE 2. HYPERBOLIC 
CROSSFLOW 
POSITIVE V 



CASE 3. HYPERBOLIC CASE 4. PARABOLIC 

CROSSFLOW CROSSFLOW 

NEGATIVE V 

Fig. 1 Eigenvalue structure in (£ 9 ) plane. 
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IB-MI = 


J(pK) #r -J(ptO* r 


-X 

0 


0 

-X 


= 0 


Solving for X from Hq. (6), one gets 


i (pi/)* + (pn* r ) ±vi (pi/)* +(pio«j , -4(pi/)* f (pn*. 

Xw ” 2(pU), 


( 6 ) 


(7) 


When Eq. (5) is satisfied, the discriminant in Eq. (7) is always 
positive. 

Now, analyzing X, and \ 2 , the following combinations are 
possible. 

1) is positive and \ 2 is negative, or X, is negative and \ 2 
is positive. 

2) \, and X^ are both positive. 

3) X/ and \ 2 are both negative. 

4) X/ or X^ is zero. 

These possible combinations are schematically shown in 
Fig. 1. Each one of these combinations describes a different 
feature of the flow in the crossflow direction tj. Referring to 
the diagrams of Fig. 1 , the following descriptions are made. 

Case 1 

Here, one eigenvalue is positive and one negative. This 
implies an elliptic-type crossflow because the characteristic 
signals are brought into point P from both the positive and 
negative direction of rj. 

Case 2 

Here, both the characteristics are positive, which means the 
characteristic signals propagate into point P only from below 
and anything happening above point P does not influence that 
point. This describes a hyperbolic-type crossflow point with a 
positive contravariant velocity K. 

Case 3 

Here, the characteristic signals propagate from above into 
point P and, similar to case 2, this describes a hyperbolic-type 
crossflow with a negative contravariant velocity V. 

Case 4 

Here, one of the eigenvalues is zero and describes a 
parabolic-type crossflow. This will represent the crossflow 
sonic line. 

The transition from an elliptic to a hyperbolic crossflow 
type takes place through a parabolic point, which is indicated 
by one of the eigenvalues going to zero. Thus, by monitoring 
the eigenvalues X, and X^, one can precisely model the 
crossflow plane terms. Depending on whether it is elliptic or 
hyperbolic, appropriate finite-difference models for the tj- 
derivative term in Eq. (1) are chosen. This will be described 
later in this paper. 

One can readily see from Eq. (7) that one of the eigenvalues 
goes to zero when =0. From the definition of p and V 

from Eq. (2), one can write 


similarly, \ 3 and \ 4 of A~ { C can also be used to determine 
the marching step size Af from a given Courant number. 


Af=min 


t 


CFL*Ai? 


CFL*An 


(9) 


The quantities (X^)^ and (X^)^ define the maximum of 
(\,,\ 2 ) and (X Jt X ¥ ), respectively, and CFL is the user- 
prescribed Courant number, usually set to values much 
greater than one for implicit schemes (CFL - 5-20). 


B. Treatment of (pU/J) { in Eq. (1) 

The direction f has been identified to be the hyperbolic 
marching direction satisfying the condition given by Eq. (5). 
Referring to Fig. 2, this derivative term will be backward 
differenced as 







where 


( 10 ) 


&/ = ur/)' 

9 - 0 first-order accurate 
= 1 second-order accurate 

Given the velocity potential <t> information at all previous 
planes i, /—I, /- 2,..., the problem is to compute 0 at the 
current plane *+!. Equation (10) involves both the density 
and contravariant velocity at the (/+ 1) plane, and both are 
functions of 0 [Eq. (1)). In order to write Eq. (10) in terms of 
0 will require only a local linearization procedure. This is 
done as follows: 


(pC/),w*<P</)i+ [(*</)* 1 40+... (11) 

where 


(pU) 4 =P 4 V+pU 4 and A0 = 0 /+/ -0, 

Substituting for p 4 and U 4 into Eq. (11), and grouping 
various term. 


//a \( t/*\a(A0) / VV\ 


3 (A0) 

drj 


— ^7 ) (8) 

Thus, when a 22 = V 2 /a 2 occurs, the method will anticipate a 
switch in the character of the crossflow and realize the onset 
of the formation of a supercritical crossflow. 

Besides providing valuable information regarding the type 
of crossflow, the eigenvalues X, and X^ of A~ l B and. 


The above locally linearized equation involves only A0 as the 
unknown to be solved for. To maintain the conservative 
differencing, both (pt/), +/ and (pt/) t appearing in the first 
term of Eq. (10) will be linearized. That is, (pU) , will be 
linearized about (/—l) plane values. The upwind differencing 
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of the {“-derivative term as shown in Eq. (10) will produce a 
truncation term whose leading term is ( 1 ~{a 2 a n /U 2 )\ 
t/^^Af. This will always represent a positive artificial 
viscosity as long as the marching condition dictated by the 
characteristic theory, Eq. (5), is satisfied. 

C. Treatment of the Crossflow Term (p Vjj) n in Eq. (I) 

Similar to the treatment of the (pU) term, the (pK) l+/ will 
also be linearized as 


(pn/ + / = (pK) / +[(pP)J l A^ + ... 

17 UV\ dA<f> / V 2 \ 3A0 

),!F + 

/ yyy\ saa „i 
+ +v ‘] 


(13) 


The above linearized expression for {pV) U! will be plugged 
inside the ij-derivative term of Eq. (1). It involves only A <t> as 
the unknown variable. The finite-difference model for the 
Ip ( y/J) ] , term will be dictated by the theory of characteristic 
signal propagation as described in Sec. II. A. When the 
eigenvalues of A ~ l B represent case l in Fig. 1 (one positive 
and one negative eigenvalue representing an elliptic type), 
then all of the terms in [p( V/J) } ? will be central-differenced. 
For this case, [a 22 - (V 2 /^)] is positive, and central dif- 
ferencing of the [p{VJ J)] f term along with the backward 
differencing of the [p(f/A/) ] f term as in Eq. (10) will 
preserve the diagonal dominance. For cases 2 and 3 of Fig. 1, 
the crossflow behaves like a hyperbolic type, and [a 22 
-{V 2 /^)] is negative. Then, central differencing of the 
terms in Eq. (13) is inappropriate, as it will destroy the 
diagonal dominance, and, in addition, will not provide the 
necessary artificial viscosity to avoid the formation of ex- 
pansion shocks. Thus, when X, and \ 2 are both positive or 
both negative (hyperbolic type), the terms in [p(V/J) ] 
should be upwind differenced depending on the direction of 
V. However, such an upwind differencing in the i) direction 
will not give rise to a tridiagonal system and, in general, the 
overall system will be pentadiagonal in nature. In order to 
preserve the tridiagonal nature of the implicit scheme, rather 
than upwind differencing the <t> derivatives, the density biasing 
concept 7,12 is implemented when the crossflow is hyperbolic. 

The procedure is as follows: 


(py) — ^ [y ^‘ a 22^ > n + ) ] 0*) 

Here, the density p has been replaced by p defined to be 
(referring to Fig. 2) 

Pi+/j+ «.* = ( ^ “ y iJ+ «.* n.k 

+ Pj- 1 +2m,k ) 0 ^) 

where m ~ 0 when ^.+^*>0,= + 1 when V u +n tk < 0. The 
artificial viscosity coefficient v f j+ ^ k is computed as follows: 


N V ' iJ+H.k 
where a is the local speed of sound and 

ii — 0 for (a„ - ^7 ^ >0 (elliptic crossflow) 

V a*'u*n.k 

= 1 for (a 22 - -7) <0 (hyperbolic crossflow) 

' a ' IJ+ H,k 


(16) 


derivative terms are central differenced in Eq. (14). Treatment 
of the density as represented by Eqs. (15) and (16) would 
always produce a positive artificial viscosity when the 
crossflow is hyperbolic. The local total velocity is always 
assumed to be greater than the speed of sound, otherwise the 
marching procedure would fail. 

In Eq. (15), the evaluation of p* depends on whether the 
flow is conical or nonconical. For conical flows, all p * 
quantities are evaluated at the /th plane. For nonconical 
flows, at each nonconical marching plane, initially p* is set to 
be the value at the /th plane and then subsequently iterated to 
convergence by setting p* to the previous iterated value of p at 
the current /+ 1 plane. 

A similar density biasing procedure is implemented for the 
[p( WfJ) ] { term in Eq. (1). 

Activating the density biasing based on the eigenvalue 
structure of A ~ l B and A ~ l C has proven to be very efficient 
in predicting sharp shock profiles. The same concept can also 
be employed for transonic applications. 


D. Implicit Factorization Algorithm 

Combining the various terms of Eq. (1) as represented by 
Eqs. (10, 14, and 15) together with the terms arising from 
[p( WlJ) ] t will result in a fully implicit model. This is solved 
using an approximate factorization implicit scheme. After 
some rearrangement of the terms, the factored implicit 
scheme becomes 

r f , Aj_ ± , / ± (p a s, \ 1 ^ p<*33 j 1 

L MBS 0 at \J At j 0dS J d£ J 


r A 2 d 1 d 

x /+ ~~ — + - — 

L /3Af dt j 0 dr) 



+ 


L d d 1 

0 dr) J 3t; J 


A$ = R 


(17) 

which has the form 

L £ L,(A0)=/? (18) 

and is implemented as 


L f ( A0)*=/?, L n (A<p) = (A<t>)\ « = + A0 (19) 

The various quantities appearing in Eq. (17) are given by 



and the right-hand-side term R consists of various known 
quantities. The algorithm Eq. (19) requires only scalar 
tridiagonal inversions. 


III. Grid System 

The transformation of the physical space (x,y,z) to a body- 
fitted computational space (f,jj,$) is performed numerically 
by using the elliptic grid generation technique of Ref. 11. The 
body geometry at every marching plane is prescribed along 
with a suitable outer boundary where freestream conditions 
are imposed. Since the equation is cast in a general coordinate 
system, the marching plane (constant f) can either be a 
constant x plane or a spherical (constant r) plane as long as the 
marching criterion (Eq. (5)) is satisfied. Given the geometry 
shape and the prescribed outer boundary, the following set of 
elliptic equations are solved to generate the interior grid. 

+ 7,, + i7 u = Q($,T7) (21) 


Thus, the density biasing is switched off smoothly when the The forcing terms P and Q are properly chosen to achieve two 

eigenvalues and X^ exhibit an elliptic crossflow. All the <£- main desirable features: 1) to cluster grid points to a bound- 
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ary, and 2) to force grid lines to intersect the boundary in a 
nearly orthogonal fashion. 

Once the grid is generated, all the metric terms a u in Eq. (2) 
and the Jacobian J in Eq. (3) are computed by numerical 
differentiation. To subtract out any numerical truncation 
error in the freestream due to incomplete metric can- 
cellation, 13 it is essential to add the term (especially for a high 
stretched nonorthogonal grid) 


d /P.t/.' 

^ a /p.k.' 

W a ( p - w ~) 

ar v / - 

s 

i 

I* 

<T5 

h 

' + V J / 


( 22 ) 


to the right-hand side of the finite-differenced model of Eq. 
(1). To be consistent with the implicit operator, Eq. (18), the 
linearization procedures given by Eqs. (12) and (14) are also 
applied in evaluating Eq. (22). 


IV. Results 

A series of calculations were performed for conical and 
nonconical geometries at various Mach numbers ( M -2-6) 
and angles of attack (a -0-10 deg) to validate the full- 
potential characteristic switch methodology and assess the 
feasibility of using numerical grid solvers for complex con- 
figurations. The results from this study are compared with 
experimental data and Euler simulation. 

The generality of the formulation allows one to choose any 
f as the marching direction, provided the condition given by 
Eq. (5) is satisfied. Thus, depending on the geometry 
definition and the flowfield character, one could choose either 
a constant jc-plane marching or constant r-plane spherical 
marching. 



ON.OFF DENSITY 

BIASING ACTIVATOR 
BASED ON 
CHARACTERISTICS 




Fig. 3 Effect of density biasing activator on the crossflow Mach 
number distribution in the shock region. 


The effect of on/off density biasing based on characteristic 
signals in the crossflow plane (ij.O, described by Eq. (16), is 
demonstrated in terms of the crossflow Mach number ( M CF ) 
distribution in the shock region in Fig. 3. When density 



Fig. 4 Grid arrangement in the marching plane for a conically 
cambered wing-body combination. 
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Fig. 5 Surface pressure distribution on a conically cambered wing- 
body combination, M ^ = 2. 
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Fig. 6 Surface pressure distribution on a flat conical wing-body 
combination, » 2. 
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biasing is applied everywhere, 7 including elliptic crossflow 
points, it introduces unnecessary artificial viscosity and tends 
to smear the discontinuities like shocks in the flowfield. This 
is seen by the dashed-line crossflow Mach number distribution 
across the bow shock and across the embedded shock on a 
cone surface in Fig. 3. When the density biasing is switched 
off at crossflow elliptic points, the shocks appear as a sharper 
discontinuity (usually within two mesh intervals), as shown by 
the solid line distribution in Fig. 3. All the calculations to be 
presented here were achieved using the second-order accurate 
implicit scheme [0= 1 in Eq. (10)1, with on/off density biasing 
activator n in Eq. (16). 

Figure 4 shows the grid arrangement in the marching plane 
for a conically cambered wing-body combination. The elliptic 
grid solver with orthogonality constraints near the surface 
required 40-60 iterations to converge to within 10"* error in 
the residual. Figure 5 shows the pressure distribution at 
A/«=2 and angles of attack of 7.81 and 10.82 deg. The 
leading-edge sweep is moderate (57 deg), and spherical plane 
marching is implemented (instead of jr-plane marching) to 
avoid low supersonic Mach number components along the x 
direction near the leading edge. The results are compared with 
experimental data given in Ref. 14. The comparison is ex- 
cellent. The marching step size Afis chosen by monitoring the 





Fig. 7 Pressure distribution on Sears-Haack body at M m -6, 
windward plane of symmetry. 



_L=r 


Fig. 8 Top and side views of a typical arrow wing-body con- 
figuration. 


eigenvalues and setting the Courant number to about 20. The 
numerical formulation, being a conservative form, predicts a 
stronger crossflow recompression on the leeward side than 
those seen in experiments. On a 20x49 (*,{) grid, the method 
requires about 1 min of CDC 7600 time. The conical flowfield 
is assumed to have converged when the change in root-mean- 
square density between two successive marching planes is 
reduced to less than 1 0 ~ 3 . 

Figure 6 shows the surface pressure distribution on a fiat 
conical wing-body (that is not designed to weaken the 
crossflow shock formation) at two different angles of attack 
(1.72 and 5.71 deg) and a Mach number of 2. The ex- 
perimental data and the numerical prediction are in excellent 
agreement and clearly indicate the presence of an embedded 
crossflow shock. 

Even though the full-potential theory is restricted by the 
isentropic assumption, one will be surprised to find that the 
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Fig. 9a Grid arrangement and surface pressure distribution for a 
symmetric arrow wing at x/t = 0.3. 
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Fig. 9b Arrow wing pressure distribution at x/t= 0.5. 




theory can be effectively utilized to predict even very high 
Mach number flows as long as A/5 is less than or of the order 
of 1 (A/5 5 1). This is demonstrated in Fig. 7, which shows the 
results for a Sears-Haack body at A/ w = 6 and different angles 
of attack (0, 4, and 8 deg). The numerical prediction is 
compared with unpublished NASA-Langley data, and the 
agreement is excellent. Constant *-plane marching is im- 
plemented for this configuration. 

Figure 8 shows a schematic of a symmetric arrow wing- 
body configuration. The actual geometry shape is prescribed 
analytically as detailed in Ref. 15. A series of computer runs 




Fig. 9d Arrow wing pressure distribution at x/t = 0.8. 


were made for this configuration at different Mach numbers 
and angles of attack and some results are presented here. 
First, an initial data plane near the nose region of the con- 
figuration is established by assuming a conical nose shape. 
The nonconical marching is then initiated. At each nonconical 
marching plane, the density is iterated to convergence (p* in 
Eq. (15) usually takes 2-3 cycles to converge to 10 ~ 5 error 
tolerance] before proceeding to the next marching plane. The 
grid at each marching plane is generated using the elliptic grid 
solver. Figures 9a-d show a series of results at x/t of 0.3, 0.5, 
0.65, and 0.8, respectively. The full-potential results are 
compared with the experimental data and Euler simulation in 
Ref. 15. Figure 9a, which shows results for an x/t of 0.3, 
clearly demonstrates the accuracy of the full-potential 
simulation. It is surprising to see that the present full- 
potential method compares with the experimental data even 
better than the Euler calculation, even at a high Mach number 
of 4.63. Similar excellent full-potential results are shown in 
Fig. 9b for an x/t of 0.5 and compared with data from Ref. 
15. The striking full-potential results are shown in Fig. 9c, 
where the unphysical oscillations experienced by the Euler 
simulation at A/* =2.36 near the wing-body junction area are 
not seen in the present method, and comparison with ex- 
perimental data is more dramatic. Figure 9d shows the 
pressure distribution at an x/t of 0.8, where the wing is 
separated from the body. The wake is simulated by assuming 
a planar shape, and imposing pressure equality (in the present 
method, it will be density equality due to full-potential for- 
mulation) across the cut. Again, the full-potential results are 
in good agreement with the Euler solution and experimental 
data. 

Figure 10 shows an angle-of-attack case, M =4.63, a = 3 
deg for the same symmetric arrow wing-body configuration. 
The results are compared with the data of Ref. 15 at x/t of 
0.65, and the agreement is good even near the wing-body 
function region. 

A typical arrow wing-body calculation using a 20 x 49 grid 
in the (»;,£ ) plane and a marching step size Courant number of 
3-5 [for a given Courant number, the predicted marching step 
size from Eq. (9) will decrease with decreasing freestream 
Mach number], required approximately 2-3 min of CPU time 
for the entire calculation. This includes the numerical grid 
generation at each plane and the conical initial data plane and 
represents an enormous savings in computer execution cost 
over other nonlinear methods, especially Euler solvers. 
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V. Conclusions 

A nonlinear full-potential aerodynamic prediction 
capability based on a sound mathematical theory of 
characteristic signal propagation has been developed. The 
method uses a general body-fitted coordinate system and 
numerical mapping techniques. The on/off density biasing 
activator in the crossflow plane has proved to be very effective 
in capturing sharp shock profiles. Results for conical and 
nonconical flows at various Mach numbers and angles of 
attack are shown to be in excellent agreement with ex- 
perimental data and Euler results. The enormous savings in 
computational cost exhibited by the present approach makes 
it a very promising substitute for the less accurate linearized 
panel methods and expensive Euler solvers, for use as a 
preliminary design tool. Future work will involve automatic 
grid generation for wing-body-nacelle-canard configurations 
and better wake treatment. 


Appendix — Derivation of the Marching 
Condition [Eq.(5)] for a Cartesian System 
The Cartesian system analog of Eq. (4) is given by 

Af x +Bf y + C/j = 0 (Al) 

where 




Equation (A I) is hyperbolic with respect to the x direction if 
1) A exist, and 2) A~ l (a£+/3C) must have real eigen- 
values for all a and 0 satisfying a 2 +/3 2 = 1 . 

Since A is assumed to exist, the eigenvalues of A~ l 
(or£ + 0C) can be obtained by setting the following deter- 
minant to zero. 

\aB+0C-\A\=O (A2) 

Substituting for A, £, and C from Eq. (Al) into Eq. (A2), the 
roots of the equation are obtained. 





(A3) 


Equation (A3) will have real values as long as the square root 
term is real. This implies the quantity inside the square root 
must be positive. Simplifying the quantity inside the square 
root, the condition becomes 
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A nonlinear method based on the full potential equation in conservation form, cast in an arbitrary coordinate 
system, has been developed to treat predominantly supersonic Rows with embedded subsonic regions. This type of 
flow field occurs frequently near the fuselage/canopy junction area and wing leading-edge regions for a 
moderately swept fighter configuration. The method uses the theory of characteristics to accurately monitor the 
type-dependent flowfield. A conservative switching scheme is developed to handle the transition from the 
supersonic marching algorithm to a subsonic relaxation procedure, and vice versa. An implicit approximate 
factorization scheme is employed to solve the finite differenced equation. Results are shown for a few configura- 
tions, including a wing/ body /wake realistic fighter model having embedded subsonic regions. 


I. Introduction 

N ONLINEAR aerodynamic prediction methods based on 
the full potential equation are used regularly for treating 
transonic 1,2 and supersonic 3 ' 5 flows over realistic wing/body 
configurations. The transonic algorithms 1,2 are designed to 
treat predominantly subsonic flows with pockets of supersonic 
regions bounded by sonic lines and shocks. The supersonic 
methods 3 * 5 are based on a marching concept and require the 
flow to remain supersonic in a given marching direction. Once 
the marching direction velocity becomes subsonic, the domain 
of dependence changes and a pure marching scheme 3 * 5 will 
violate the rules of characteristic signal propagation. The 
possibility of a marching velocity becoming subsonic in a 
supersonic flow is great, especially for low supersonic 
freestream Mach number flows (A/^ = 1.3 - 1.7) over mod- 
erately swept fighter-like configurations (sweep angle A = 45 — 
50 deg) and over forebody shapes having a sizeable 
fuselage/canopy junction region. There is a strong need to 
construct a supersonic marching computer program that has 
built-in logics to detect and treat the embedded subsonic regions. 

The method of Ref. 5 is based on the characteristic theory 
of signal propagation and uses a generalized, nonorthogonal, 
curvilinear coordinate system. Compared to other nonlinear 
supersonic methods, 3 the method of Ref. 5 has no restrictions 
(limitations of the full potential theory hold) on its applicabil- 
ity to complex geometries and intricate shocked flowfields. It 
is a conservative formulation and uses numerical mapping 
techniques to generate the body-fitted system. The purpose of 
this paper is to describe an extension to the methodology of 
Ref. 5 to include the treatment of embedded subsonic regions 
in a supersonic flow. 

The paper describes the characteristic theory involved in 
determining the condition for a marching direction to exist. 
Once that condition is violated, the marching scheme is transi- 


Received June 30. 1983; presented as Paper 83-1887 at the AIAA 
Computational Fluid Dynamics Meeting, Danvers, Mass., July 13-15, 
1983; revision received March 19, 1984. Copyright © American In- 
stitute of Aeronautics and Astronautics Inc., 1984. All rights reserved. 

* Manager, Computational Fluid Dynamics Group. Associate Fel- 
low AIAA. 

fMember, Technical Staff. Member AIAA. 

^Professor, Department of Mathematics. Member AIAA. 


tioned to a relaxation scheme through a conservative switching 
operator. For marching condition violation, the total velocity 
q does not have to be subsonic. Even for a supersonic total 
velocity q , if the component in the marching direction is 
subsonic, a relaxation scheme is required. In order to properly 
produce the necessary artificial viscosity through density bias- 
ing, the paper defines two situations: 1) the total velocity q is 
supersonic, but the marching direction component is subsonic 
[defined as marching subsonic region (MSR)]; and 2) the total 
velocity q is subsonic [termed as total subsonic region (TSR)]. 

Results are presented for a few configurations that exhibit 
either the MSR or both the MSR and TSR flowfield. The 
paper also presents results from a wake model applied to a 
realistic wing/body fighter configuration. 

The Appendix describes a flux biasing concept that will 
supersede the density biasing procedures currently in use. 

The methodology of this paper is not restricted to the full 
potential equation alone. Currently, similar marching/ relaxa- 
tion methods are under development at Rockwell for applica- 
tion in parabolized Navier-Stokes (PNS) codes to treat the 
embedded subsonic regions or streamwise separated flows 
without having to use a time-dependent Navier-Stokes pro- 
gram. 


II. Equation and Characteristic Theory 

The conservative full potential equation cast in an arbitrary 
coordinate system defined by £ =» f(jc t y, z), tj = Tj(x,y, z), 
and £ = £(x,y,z), takes the form 

k vwm.- <■> 

where U, V, and W are the contravariant velocity compo- 
nents. Introducing the following notation for convenience: 



F= U 2 , 

II 

X = X/, 

II 

z ~ x i 

'■'-T 

II 

V = X 2 , 
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the contravariant velocities and density are given by 


Eq. (5) are given by solving the quadratic 


Ui-Zo.jtxj 1.2.3 

j-t 

* dx t dXj i — 1,2*3 / f 

a u = 2* ^ , 3 (transformation metncs) 

dx k ° x k 

p ~ [' - u *t + v *' + ~ 1 }J 

a = speed of sound = \/ML (2) 

The Jacobian of the transformation J is represented by 


/- 


B( x *y>z) 


f, f, !; 

Vx Vv V, 

i, «: 


( 3 ) 


Equation (1) is in terms of a general coordinate system (£, tj, £) 
and can accommodate any kind of mapping procedure, either 
analytical (conformal mapping) or numerical. The nature of 
Eq. (1) can be analyzed by studying the eigenvalue system of 
Eq. (1). Combining the irrotationality condition in the (£,77) 
and (f,{) plane and Eq. (1), one can write the following 
matrix equation: 

^/ f +y,+c/«-o (4) 

where 

L j{pU)>< !j(pUU, ■ jipuu, ’ 

A ~ 0 1 0 

0 0 1 


j(py)*, L j{pv)*, 7 (pv)*, 

-100 

000 


~j(pW)*t far)*, 7 (piv)*, 

000 

-100 




/= 




-\ 2 (pU)<, l + \[a(pV)<, ! + B(pW)<, s + a(pU) < ,, 

+q ( pu)*(\ ~{<*'(py)*, +^'(pw / )« < 

+aP[{pW)*,+(pV)<, l ]}-0 (6) 

Representing Eq. (6) in the form 


A\ : + B\ + C = 0 (7) 

the discriminant ( B 2 - 4AC) determines the character of Eq. 

( 4 ): _ 

1) If ( B : - 4AC ) remains positive for all a and ft satisfy- 
ing or + ft 2 = 1, then the eigenvalues of Eq. (4) are real and 
direction^* is hyperbolic (marching scheme is valid). 

2) If ( B 2 - 4A C) is negative , then the eigenvalues of Eq. (4) 
are complex and direction £ is elliptic (requires a relaxation 
method). 

To analyze when the eigenvalue solutions of_Eq. (6) are real 
and when complex, the discriminant ( B 2 - 4AC) is rewritten 
in the following form using Eq. (2): 


V uvy' 

( U 2 \ 

( y: \] 

a n — r - 

a U r 


L\ a - / 

l a- J 

l cr I 


+ 2 aft 

[h-f) 

(--?)- 

■(--7) 

(- 7)1 


uw\ : 

-(--7) 

(--7) 

] < 8 > 


Using the properties of a positive definite quadratic form and 
the Schwarz inequality {a^a^ > ajj), Eq. (8) can be shown to 
have the following results: 

1) (B 2 - 4AC) is positive if [a n -{U : /a 2 )\ is less than 
zero. Then the f direction is hyperbolic (the marching al- 
gorithm of Ref^S is valid). 

2) ( B 2 - 4 A C) is negative it [a,, - (U 2 /a 2 )] is greater than 
zero. Then the f direction is elliptic (requires a relaxation 
scheme). 

Physical Interpretation 

The physical interpretation of these results from the char- 
acteristic theory is illustrated in Fig. 1. Let q be the total 
velocity. The projection of q in the direction normal to the 
f = constant surface is given by 


q */? = («/ + vj + wk ) 


(U + £■■/+£;*) 

i; + {;■+{: 


7 


( 9 ) 


The subscripts in Eq. (4) denote differentiation with respect to 
that variable. 

The matrices A, 5, and C appearing in Eq. (4) can now be 
analyzed to determine the character of that equation. In 
general, the following is true: 

1) Equation (4) is elliptic in the f direction if the matrix 
A (aB + ftC) has complex eigenvalues for all combinations 
of a and ft such that or + ft 2 *= 1. 

2) Equation (4) is hyperbolic in the £ direction if A~ ! (aB + 
ft C) has real eigenvalues for all a and ft satisfying a 2 + ft 2 — 1. 

The eigenvalue structure of A ~ ! {aB + ftC) can be obtained 
by setting the determinant 

\aB 4- ftC - A/4| = 0 (assuming A exists) (5) 
Substituting for A, B , and C from Eq. (4), the eigenvalues of 


where u, u, and w are the Cartesian velocities and n the 
normal to the £ = constant plane. Figure la shows the case 
when V/yfai~i is greater than the speed of sound {[a n — 
(U 2 /a 2 )] c 0}. For this case, the characteristic cone of in- 
fluence is behind the £ = constant plane and marching along 
is valid. Figure lb illustrates the case for the q> a, but for the 
U/fc < a situation, [a n - (U 2 /a 2 )] > 0. For this case, a 
part of the characteristic cone of influence lies forward of the 
£■ = constant plane and marching along f is not possible. This 
case (Fig. lb) is termed marching subsonic region (MSR) in 
this paper. Figure lc shows the case when q < a and 
< a, [a ,, — ( U : /a ■)] > 0. This represents a pure subsonic flow 
and marching along f is not possible. This case is termed total 
subsonic region (TSR). For cases represented in Figs, lb and 
lc, a relaxation algorithm is required. 
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I - TOTAL VELOCITY 
i - NORMAL TO 

f - CONSTANT PLANE ' ^CHARACTERISTIC 

- SPEED OF SOUNO Z0NE 0F INFLUENCE 

U 2 

b) MARCHING SUBSONIC REGION (MSR). (a^ --)>0.q >a 
MARCHING ALONG f IS NOT VALID. a 



COMPLEX CHARACTERISTICS. 


Fig. I Role of characteristics in defining supersonic region, marching 
subsonic region (MSR), and total subsonic region (TSR). 


HI. Numerical Method 

Figure 2 shows the schematic of a fuselage/canopy fore- 
body geometry with an embedded MSR and TSR present in a 
supersonic flow. To solve this problem, the marching scheme 
of Ref. 5 will be used when [a u - {V 2 /a 2 )] is negative and a 
relaxation scheme when [a n - (U 2 /a 2 )] is positive. First, 
march from the nose up to the plane denoted by (A-B) in Fig. 
2, using the method of Ref. 5. Then, between (A-B) and 
(C-D), which embed the subsonic bubble (MSR and TSR), use 
a relaxation scheme and iterate until the subsonic bubble is 
fully captured. Then, resume the marching scheme from the 
plane (C-D), downstream of the body. 

The purpose of this paper is to present a conservative 
algorithm that will automatically switch from a pure marching 
scheme of Ref. 5 to a relaxation method at the onset of an 
MSR formation and revert to the marching procedure when 
the flow becomes fully supersonic again. The entire flowfield 
can be classified into three types with respect to the marching 
direction 

1) At a grid point, the marching direction is hyperbolic and 
the total velocity q is supersonic, [a n — (U 2 /a 2 )\ <0, q > a. 
This point will use the algorithm of Ref. 5. 

2) At a grid point, . the marching direction f is elliptic , 
[a n - (U 2 /a 2 )] > 0, but the total velocity q is supersonic, 
q > a (MSR). This point will be treated by a transonic oper- 
ator with a built-in density biasing based on the magnitude of 
[1 -{a 2 /q 2 )]. 

3) At a grid point, the direction f is elliptic and the total 
velocity q is subsonic , q < a (TSR). This point will be treated 
by a subsonic central differenced operator. 


Treatment of (3/3£)Ip(t//«DI in Eq. (1) 

Refer to the computational molecule in Fig. 3. 



Fig. 2 Embedded subsonic bubble in a supersonic How. 


where 


5 refers to backward differencing 
d refers to forward differencing 



In Eq. (10), the first term corresponds to the supersonic 
marching operator of Ref. 5 and the second term is the 
subsonic operator. 

The backward difference operator in Eq. (10) is represented 
by 


d ( U) 

1 • 5 

7 u 2 

\ 3 . 


Td p 7) 

i + / e/f 


fat** 

/ 

uv\ d 


uw\ 

d , 




-) 

a| A * + t/ < 

’ 




= 


(11) 


The term (<9/d£)(A<>) is backward differenced. Reference 5 
gives more details on this supersonic marching operator. 

The forward difference operator in Eq. (10) is represented 
by 

* A. 

where 

P n ,Zl = P% / - v( p%, - p'/), for U > 0 

v= max[o, 1 -(a 2 /q 2 )\ (13) 

The superscript n + / denotes the current relaxation cycle for 
a subsonic bubble calculation. 

Note that in Eq. (12) the term <f>~ is backward differenced 
such that (d/d$)(p/J)a t i<j>£ will provide the central differenc- 
ing needed for an elliptic (subsonic) point. The density biasing 
[Eq. (13)] is activated only when the total velocity q is greater 
than the speed of sound a. This will take place when a grid 
point is in the region denoted by MSR in Fig. 2. When q < a 
(the TSR in Fig. 2), the density is not biased and the genera- 
tion of artificial viscosity is turned off. The <t> derivatives in 
Eq. (13) can be rewritten in terms of £<p, just as in Eq. (11). 
Equation (10) can also be interpreted as 


d 1 u\ 

d j 

f U\ Ky d. 

3 / 

r u\ 

3?l p 7r 

“ d$\ 

t J /( + / 

ad 

> p 7/ 


( a n is + “i.a + “ia ) , . 


( 12 ) 


supersonic 


marching 

subsonic 


elliptic 

operator 


flux biasing to produce 
the artificial viscosity 
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Fig. 3 Conservative type-dependent switching scheme for the treat- 
ment of subsonic bubble in a supersonic flow. 


Figure 3 illustrates various possibilities that can be handled 
by Eq (10). It has both the shock point operator and the sonic 
operator required to treat the type- dependent flow. The only 
issue that philosophically affects the concept of a conservative 
scheme is that the definition of pU for a supersonic operator 
in Eq. (11) is different from the definition for the subsonic 
operator of Eq. (12). 

The evaluation of the subsonic operator in Eq. (12) requires 
velocity potential <f> values at / + 1 and i + 2 planes from the 
previous n relaxation cycle to compute the density. The 
section on initial and boundary conditions below prescribes a 
method to start the first relaxation cycle of the subsonic 
bubble calculation. 

Treatment of (d/3Tj)|p(F//)| in Eq. (1) 

Referring to the Fig. 3 a molecule. 


In Eq. (16), the evaluation of p* depends on whether the 
flow is conical or nonconical. For conical flows, all p* quanti- 
ties are evaluated at the jth plane. For nonconical flows, at 
each nonconical marching plane, initially p* is set to be the 
value at the / th plane and then subsequently iterated to 
convergence by setting p* to the previous iterated value of p 
at the current i + 1 plane. Reference 5 provides more details 
on the density biasing procedure and the implicit treatment of 
&/dimV/J)]j+L in Eq. (15). 

When the point is elliptic , the density biasing is defined by 


+ a?) 


where v — max[0,l - {a 2 /q 2 )\. As before, the superscript n +■ 

1 denotes the current relaxation cycle for a subsonic bubble 
calculation. Note the difference in the definition of v and v. 
The density biasing in the cross-flow direction 17 is turned off 
when the total velocity q is less than the speed of sound u, 
just as in the marching f direction (Eq. (13)]. The implicit 
treatment of V in the marching subsonic operator of Eq. (15) 
is the same as that of the supersonic part, explained in Ref. 5. 

A similar procedure is implemented for the [p( W/J)\ ( term 
in Eq. (1). 

Implicit Factorization Algorithm 

Combining the various terms of Eq. (1) as represented by 
Eqs. (10)-(17) together with the terms arising from [p(1F/J)] € 
will result in a fully implicit model. This is solved using an 
approximate factorization implicit scheme. After some re- 
arrangement of the terms, the factored implicit scheme be- 
comes 


d f d 1 d ( p a u\ l d P a u d 

JIi;Jl + pTz\Ju) + i}Ti~rTz 


a, a 

,L±( 

' P°:i \ 

1 d pa„ 


m 3v H 

P 3n\ 


+ P dr, J 

dj) 


(18) 

The density p appearing in Eq. (18) can be either p or p 
depending on the sign of [a n - (U : /a 2 ) ] as illustrated in Eq. 
(15). 

Equation (18) has the form 




+(/- 



supersonic 


marching 

subsonic 


05) 


L ( L,(A<t,)-R 

and it is implemented as follows: 
L ( (A<t>)' = R L v ( A*) ■*(&<(,)' 


(19) 

<f> t y = </>, + A<£ 

( 20 ) 


where 



if ! 


<0 

1 + / 

(supersonic point) 

*,./ = 0 

if | 

h-5] 

1 >0 

1 + / 

(MSR) 


The various quantities appearing in Eq. (18) are given by 



When 0 I + J = 1, that is, the point is supersonic with respect 
to f, only the first term in Eq. (15) is used and the biased 
density p is defined by (for V > 0), 

+P%l) ( 16 ) 



where 


+ „ = ( 21 ) 


► “ max 


1 0,1 -a 2: 



and the right-hand side term R consists of various known 
quantities. 
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If the flowfield does not contain an embedded MSR or 
TSR, the implicit factored algorithm of Eq. (18) performs a 
pure marching procedure starting from an initial known data 
plane. In this situation, there is no need to go back to the 
upstream starting plane and iterate the solution. However, if a 
subsonic bubble is present (between planes AB and CD in 
Fig. 2), then the solution procedure of Eq. (18) performs a 
relaxation method and iterates for the elliptic subsonic bubble 
to converge (superscript n in Eqs. (12), (13), and (17) refers to 
the relaxation cycle counter). 

Initial and Boundary Conditions 

Initial Conditions 

For a pure supersonic flow, initial conditions need to be 
prescribed only at the starting plane. Usually, the starting 
plane is set close to the apex of the configuration to be solved 
and the conical solutions are prescribed. 

Inside an MSR, as in Fig. 2, when Eq. (12) is applied at an 
(i + 1) grid point, information on <f> i+2 is required to form the 
density p and various derivative terms. For the first relaxation 
pass, an initial estimate for quantities in the (i + 2) plane is 
prescribed in the following manner: 


JK 



MSR operator 


needs initial 
estimate 


- ( j ( a " *r + “'A + a / A ))} ( + 


( 22 ) 


In Eq. (22), sonic conditions are assumed at (i + 2) for the 
first relaxation pass. 


^ + :“<7*V f ( *//), + > (23) 

The sonic values p* and q* are purely a function of the 
freestream Mach number M x . Also, p, + / in Eq. (22) is 
initialized to be p,. 

For the second relaxation cycle and onward («>1), the 
conditions from the previous relaxation cycle are used, 



(24) 


Boundary Conditions 

At a solid boundary, the contravariant velocity V is set to 
zero. Exact implementation of V - 0 in the implicit treatment 
of Eq. (18) is described in Ref. 4. 

The outer boundary is set away from the bow shock and the 
freestream velocity potential is imposed along that 
boundary. All discontinuities in the flowfield are captured. 
The precise, density biasing activator y, based on the char- 
acteristic theory, allows for sharp capturing of shocks in the 
flow. 

Behind the trailing edge of a wing, a wake model is im- 
posed. Figure 4 shows a schematic of a wake model. At a 
point P lying on the wake, the boundary condition is that 
there is no jump in the pressure across the wake, i.e., ( p p - 
Pq ) = 0. In the full potential (isentropic) formulation, this 
translates into the condition that the jump in density (p p - p Q ) 
is zero, or the jump in the total velocity q is zero [(q P - q Q ) = 
0]. The jump in q across the wake is set to zero in an 
approximate manner in the following way. 

First, compute the jump in the potential <f> at the trailing- 
edge point P* and maintain that jump constant along the line 
P'P in Fig. 4. At the wake point P, Eq. (1) is not valid. 
Instead of solving Eq. (1), <p nv = 0 is satisfied at the wake 
point P to achieve the condition (^)/> - (<>,)^ = 0. Incorpo- 
rating a constant jump in along P'P insures (<£ f )/> - (fy ) Q 



, Vp- , % =0 * ( °-M>V 0 

Fig. 4 Wake boundary condition. 



Fig. 5 Axial surface pressure distribution for a developed cross-sec- 
tion forebody a 1.7, a” -5 deg). 

=?0. The net effect is that (q P ~ q Q ) is approximately set to 
zero, yielding the necessary wake boundary condition. The 
following section presents a calculation performed for a realis- 
tic wing/body/wake fighter model and shows an excellent 
matching of the pressures across the wake, using the above 
wake boundary condition. 

Grid System 

The transformation from the physical space (*,>\z) to a 
body-fitted computational space (f, tj, £) is performed numeri- 
cally at each constant { plane by using the elliptic grid 
generation technique of Ref. 6. Once the grid is generated, all 
the metric terms a tJ in Eq. (2) and the Jacobian J in Eq. (3) 
are computed by numerical differentiation. As described in 
Ref. 5, a freestream error subtraction is performed at each grid 
point to account for any improper metric cancellation. 

Density Biasing Summary 

This section summarizes in a tabular form the type-depen- 
dent density biasing procedures incorporated in this paper to 
generate the proper artificial viscosity. See Table 1. 
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Table 1 Summary to type-dependent density biasing procedure 


Term 

Total supersonic 

Marching subsonic 

Total subsonic 

Definition 

(a,, -({/-’/«- »<o, 

q> a 

[a u -{U 2 /a : ) ] > 0. 
q > a 

[a n ~{U 2 /a')\ > 0, 
q< a 

pU in 

f direction 

Upwind differencing, 

Eq. (ID 

Density biasing 
based on 

[\-{a 2 /q 2 )\ in Eq. (13) 

Shut off density 
biasing 

pV,pW in 7j,£ 
directions 

Density biasing 
based on 

[1 -a 22 (a 2 /V 2 )). 
[\-a )3 {a 2 /W 2 )\ 
p in Eq. (16) 

Density biasing 
based on 
[1 -< a-/*')] 
p in Eq. (17) 

Shut off density 
biasing 



Fig. 6 Circumferential pressure distribution for a developed cross-sec- 
tion forebody (Af*. = 1.7, a= -5 deg, x// = 0.28). 


IV. Results 

As illustrated in Fig. 2, supersonic marching calculations 
are performed from the nose until an embedded MSR forms. 
In Fig. 2, the plane AB is the last supersonic marching plane 
preceding the subsonic bubble and forms the upstream com- 
putational boundary for the relaxation calculation. For the 
first relaxation pass through the subsonic bubble region, 0, + / 
in Eq. (10) is set equal to 0, and (pt/)/ + » “ P*< 7 *. From the 
second relaxation cycle on, 0 I + / , $ n and ( f>U) i+2 are com- 
puted according to their definitions. A typical supersonic flow 
with a subsonic bubble calculation required at most only four 
relaxation cycles (iterating back and forth between planes AB 
and CD) to obtain a converged location’ for the bubble. The 
initial guess, based on the sonic conditions p*q *, worked out 
very well for all the subsonic bubble cases presented in this 
paper. The ( 77 , £) marching plane can be any arbitrary surface, 
but for convenience was chosen to be a constant x plane. 

The step size in the marching direction f for the supersonic 
part {[a n - (U 2 /a 2 )\< 0 , q> a) was automatically chosen 
by setting the Courant number 5 to be around 5. Once the 
MSR forms, the eigenvalues become complex and the step size 
cannot be computed based on a specified Courant number. 
For marching planes containing the MSR/TSR, the step size 
was specified into the code depending on the geometry varia- 
tion. When geometry changes were drastic (region of emer- 
gence of a wing from a fuselage), usually a smaller step size Af 
was required (as small as 0.003 - 0.005 for a total length of 
one) to properly account for rapid changes in the flow. Once 
the MSR/TSR is fully captured and the flow becomes super- 
sonic again, the step size selection once again becomes based 
on the Courant number. For a pure supersonic flow all the 
way, the entire calculation could be performed using 40 planes 



*(in) 

Fig. 7 Nose region geometry for Space Shuttle. 



Fig. 8 Surface pressure distribution at leeward plane of symmetry. 


or less (A£> 0.025). However, once an MSR or TSR is 
present, the total number of £ planes in the calculation could 
go as high as 300. 

Figure 5 shows the surface pressure distribution in the axial 
direction on the upper (0 = 0 , lee side) and lower (0 = 180 
deg, windward side) plane of symmetry for a developed cross- 
section forebody geometry reported in Ref. 7, At A/ x = 1.7 
and a = - 5 deg, the lee side has an embedded MSR that 
required use of the relaxation operator in Eq. (10). A pure 
supersonic x marching for this case would have failed without 
the MSR treatment described in this paper. 

Figure 6 shows the circumferential pressure distribution for 
the same developed cross-section forebody at A/ x = 1.70, a = 
“5 deg, and x// = 0.28. The embedded MSR thickness is the 
largest at this axial station. The extent of the subsonic bubble 
is marked in Fig. 6 . The results of Figs. 5 and 6 exhibit only 
MSR — TSR is not present. 

To simulate both the MSR and TSR, the flow over the 
Shuttle orbiter at A/ x = 1.4 and a = 0 deg was considered. 
The side view, cross section, and grid in the fuselage/canopy 
region of the orbiter are shown in Fig. 7. At this Mach 
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Fig. 9 Supersonic fighter with an embedded marching subsonic region 
near the leading edge. 


C 


P 


*“ * * 0.45 

Fig. 10 Pressure distribution on a fighter-like configuration ( W x “ 1.6, 
a ** 5 deg). 



number, the fuselage/canopy junction exhibits a large MSR/ 
TSR. Figure 8 shows the surface pressure distribution along 
the leeward plane of symmetry. At x s 4.3 m (170 in.), which 
is the beginning of the canopy, the pressure increases rapidly 
from C p s 0.3 to 1.0 and an MSR/TSR is formed. It required 
three relaxation cycles to develop the solution. The compari- 
son with the Rockwell experimental data is favorable. The 
blunt body initial solution for this Shuttle case was obtained 
from the unsteady full potential code of Ref. 8. 

Figure 9 shows a supersonic fighter configuration with a 
wing sweep of around 48 deg. At a freestream Mach number 
of 1.6 and a - 5 deg, the leading edge of the wing exhibits an 
MSR/TSR. To solve the flowfield over such a fighter con- 



Fig. 11 Grid and pressure distribution in the wake region of a 
fighter-like configuration (Af^ = 1.6, a « 5 deg, jr = 0.85). 



Fig. 12 Circumferential pressure distribution in the vertical tail and 
wing region of a fighter-like configuration = 1.6, a = 4.46 deg, 
x // =* 0.90). 


Table 2 Test cases for fighter-like configurations 


a, deg 

5 

5 

5 

5 


1.6“ 

1.6 b 

1.4 b 

1.6 b 

A, deg 

48 

48 

48 

55 

Q 

Code 

0.298 

0.3016 

0.3561 

0.29186 

Data c 

0.277 

0.295 

0.342 

0.3 

Q> 

Code 

0.0462 

0.04916 

0.04117 

0.028129 

Data* 

0.0457 

0.0493 

0.0425 

0.0301 

•Tail off. 

b Tail on. 

'Rockwell data. 




figuration, one needs to use the embedded subsonic bubble 
treatment. Figure 10 shows the surface pressure at various 
axial stations along with respective grid distribution for the 
wing/body geometry. For this case, the MSR/TSR starts 
around x = 0.4. Figure 11 shows the pressure distribution for 
the fighter configuration of Fig. 9 at an axial station x/l — 0.85, 
where a wake sheet is present. The grid distribution goes 
around the wake sheet just like a wing/body case. The ap- 
proximate wake model described in the paper seems to pro- 
vide the correct zero pressure jump condition across the wake, 
as seen in Fig. 11. Figure 11 shows the simulation without the 
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Consider the term {d/dT))[p{V/J)\ in Eq. (15). The density 
biasing procedure defines p to be 

f,* ! = ( / - >0 p, * t + M p, + Pj- 1) (ai) 

where 

v = max[ 0, 1 ~{a 2 /q 2 )\ 


In the flux biasing technique, it will be modified to 


(A2) 


where 


Fig. 13 Drag prediction for a double-wedge delta wing at =» 1.62, (pq)_=*0 it q < a 

a « 0 deg; for sweep angles less than 60 deg, the leading edge has a ' 

marching subsonic flow (C^ - C D/ + Cn L UU )- = (pq) ~ P*q* if q> a. 


vertical tails. Figure 12 shows the result for a different fighter 
model with a pronounced wing/body shape and a vertical tail. 
At this cross section, x//=»0.9, the geometry is multiply 
connected with a wake sheet present between the tail and the 
wing. The circumferential pressure distribution on the wing/ 
body/tail/wake and the gridding are shown in Fig. 12. 

The lift and drag coefficients from the present calculation 
for a fighter model are given in Table 2. The comparison with 
Rockwell experimental data is excellent. 

Figure 13 shows the drag prediction capability of the full 
potential code by demonstrating it on a double- wedge delta 
wing at =* 1.62. At this Mach number, the leading edge 
exhibited the presence of an MSR for sweep angles less than 
60 deg. A pure supersonic marching code would not have 
worked for this case. The drag calculation from the full 
potential code compared very well with the experimental data 
available in the Princeton series. 

V. Conclusions 

A nonlinear full potential method has been developed to 
treat supersonic flows with embedded subsonic regions. A 
conservative switching scheme is employed to transition from 
the supersonic marching algorithm to a subsonic relaxation 
procedure. The theory of characteristic signal propagation 
plays a key role in activating various density biasing proce- 
dures to produce the necessary artificial viscosity. The method 
has been shown to produce results that were hitherto not 
possible using a pure supersonic marching scheme. The con- 
cept of density biasing will be modified in the future to a flux 
biasing procedure described in the Appendix. 

Appendix: Flux Biasing Procedure 

Based on the work of Hafez et al., 9 it is possible to modify 
the density biasing concept to a flux biasing procedure. 


where p* and q* represent sonic conditions, a the local speed 
of sound, and (pq) the flux. When the flow is purely subsonic, 
the flux biasing is turned off automatically. 
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An aerodynamic prediction technique based on the steady form of the full-potential equation has been applied 
to a variety of three-dimensional supersonic flow problems exhibiting embedded subsonic regions. A conser- 
vative switching scheme Is employed to transition from the supersonic marching procedure to a subsonic relaxa- 
tion algorithm, and vice versa. Numerical solutions are obtained for a number of complex configurations, in- 
cluding advanced tactical fighter, Langley canard-wing fighter configuration, isolated shuttle orbiter, and mated 
shuttle orbiter configuration with external tank. The computed results are In good agreement with available ex- 
perimental data. 


Nomenclature 

a = speed of sound 

a,, = transformation metric + $ 

C D = drag coefficient 

C L = lift coefficient 

Cm = pitch-moment coefficient 

C p = pressure coefficient 

i t j,k = streamwise, radial, and circumferential indices 

J - Jacobian of transformation 

= freestream Mach number 
q+ =[p*>~ 1 /mj] ,/i , sonic condition 

U t V t W - contravariant velocities 

x ,y t z Cartesian coordinates 

a = angle of attack 

7 = ratio of specific heats 

= transformed coordinates 
p = density 

p* = sonic density 

^ = velocity potential 

( ) = wing sweep angle 


Introduction 

T HE prediction of inviscid low supersonic Mach number 
flowfields about complex three-dimensional configura- 
tions is of great interest to both researchers and designers. For 
treatment of such flows, full-potential methods 1 * 3 based on a 
space-marching procedure offer the advantage of requiring 
only moderate computer resources (memory and time) while 
maintaining sufficient accuracy. 

In the full-potential method of Refs. 1 and 2, the equation is 
transformed to a generalized, nonorthogonal, curvilinear 
coodinate system and is solved by a highly efficient, implicit, 
finite difference scheme based on the characteristic theory of 
signal propagation. A space-marching technique is used when 
the flow is supersonic in a given marching direction. If the 
velocity in the marching direction becomes subsonic, the do- 
main of dependence changes and the marching scheme is 
modified to a relaxation-type method through a conservative 
switching operator. 


Presented as Paper 85-0272 at the AIAA 23rd Aerospace Sciences 
Meeting, Reno, NV, Jan. 14-17, 1985; received March 14, 1985; revi- 
sion received Sept. 10, 1985. Copyright © American Institute of 
Aeronautics and Astronautics, Inc., 1985. All rights reserved. 
•Member Technical Staff. Member AIAA. 
tManager, Computational Fluid Dynamics Department. Associate 
Fellow AIAA. 
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The presence of subsonic pockets in a supersonic flow oc- 
curs very frequently near fuselage-canopy junction areas and 
wing leading-edge regions. In fact, future design of advanced 
fighter wings (A/* = 1. 2-2.0, wing sweep - 48 deg) will pur- 
posely incorporate subsonic regions near the leading edge to 
benefit from the leading-edge suction peak associated with 
subsonic flows. 

In Ref. 1, a numerical mapping technique is used to 
generate the body-fitted coordinate system at a marching 
plane. The key advantage of this method is that it has no 
restrictions on its applicability to complex geometries and in- 
tricate shocked flowfields. In contrast to the general coor- 
dinate formulation of Ref. 1, the method of Ref. 3 is based on 
a spherical plane marching technique and its application to 
general three-dimensional geometries is yet to be 
demonstrated. 

The main purpose of this study is to investigate the 
usefulness of the methodology of Ref. 1 in treating supersonic 
flows with large embedded subsonic regions over complex 
geometries, including realistic fighter configurations, shuttle 
orbiter, and multibody configurations (orbiter on top of the 
external tank/solid rocket boosters) at low supersonic Mach 
numbers (A/* - 1.2 to 2.0). 

Analysis 

The physical and computational coordinate systems are 
shown in Fig. 1. As discussed in Refs. 1 and 2, the entire 
flowfield is divided into three regions (see Fig. 2); 1) the pure 
supersonic region, 2) the marching subsonic region (MSR), 
and 3) the total subsonic region (TSR). The basic governing 
equations and boundary conditions are essentially the same as 
in Ref. 2 and, therefore, only a brief discussion of the method 
is presented here. 

Governing Equation 

The conservative form of the full-potential equation cast in 
an arbitrary coordinate system defined by f=f(x,y,z), 
l = l(x,y,z), and £ = £(x,y,z) can be written as 



where the density is given by 

P= [l [w» r + y<t>, + W<t>( - 1]] (2) 

and is the freestream Mach number, U t V t and W are the 
contravariant velocity components, and J is the Jacobian of 



1080 


SZHMA, RIBA, SHANKAR, AND GORSKI 


J. AIRCRAFT 




Fig. 2 Embedded subsonic bubble in a supersonic flow. 



Fig. 3 Fighter-like configuration (ATF). 


the transformation. The treatment of each term in Eq. (1), in- 
cluding the density biasing procedure and the implicit approx- 
imate factorization algorithm, can be found in Refs. I and 2. 

Initial Conditions 

Supersonic Flow Region 

For a pure supersonic flow, initial conditions are required at 
the starting plane. For sharp-nosed configurations, conical 
solutions are prescribed, and for a blunt-nosed configuration, 
the axisymmetric, unsteady, full-potential solver of Ref. 4 is 
used to obtain the detached bow shock flowfield in the 
forebody region. 

Embedded Subsonic Flow Region 

When Eq. (1) is applied at the (i + 1) plane within an embed- 
ded subsonic region, information on the flux pU at the (/+2) 
plane is required. For the first relaxation pass, sonic condi- 
tions are assumed at (/ + 2) 

p ( .2=p , =(-4r+ 2 ^-Afi) 1 ^-' «) 

\ 7 + 1 7+1 / 

Ui. 2=<?*'/0|l j + 2 

where 

<7* = VA/J.P 

Sonic- values p • and q* are purely a function of the 
freestream Mach number The quantity a n is a transfor- 
mation metric term. 


Table 1 Teat cases 

Case 1 Case 2 Case 3 Case 4 

( Fi g- 3 ) (Fig. 10) (Fig. 12) (Fig. 17) 

My, 1.6, 1.4 2.0 1.4 1.4 

a See Table 2 4.0 0 deg, -1.94 deg 0 deg 


Boundary Conditions 

In order to solve the full-potential equation, it is essential to 
specify appropriate boundary conditions on the body surface 
and along the outer boundary. 

Body Surface 

At a solid boundary, the contravariant velocity V is set to 
zero. Exact implementation of V=0 in the implicit treatment 
of Eq. (2) is described in Ref. 1. 

Outer Boundary 

The outer boundary is outside the bow shock where the free- 
stream velocity potential is imposed. All discontinuities in 
the flowfield are captured. The precise density biasing ac- 
tivator of 'Ref. 1, based on the characteristic theory, allows for 
sharp capturing of shocks in the flow. 

Swept Trailing-Edge Wake Treatment 

In order to treat the region behind the trailing edge, an ar- 
tificial cut is created, and the pressure jump [p] across this cut 
is imposed to be zero as a boundary condition. The full- 
potential equation is not solved at grid points on the wake cut. 
Instead, 0„=O is solved to provide [p]=0 across the wake 
cut. A complete discussion of this is given in Ref. 2. 

Method of Solution 

Figure 2 shows the schematic of a fuselage-canopy forebody 
geometry with an embedded MSR and TSR present in a super- 
sonic flow. To solve this problem, the marching scheme of 
Ref. 1 is used when (a n - Lft/a 2 ) is negative, and a relaxation 
scheme is used when (o u - LP/a 2 ) is positive. 

First, march from the nose up to the plane denoted by A-B 
in Fig. 2, using the method of Ref. 1. Then, between planes 
A-B and C-D, which embed the subsonic bubble (MSR and 
TSR), use a relaxation scheme and iterate until the subsonic 
bubble is fully captured. Finally, resume the marching scheme 
from the plane C-D downstream of the body. 

Geometry and Grid System 

The geometry of a configuration is prescribed at discrete 
points in a cross plane (usually x= constant planes) at various 
axial locations. These geometry input points are usually ob- 
tained from a geometry package such as GEMPAK 3 or CDS. 6 
The input points are then divided into several patches, and at 
each patch a key-point system is established. The geometry at 
a marching plane is then obtained by joining appropriate key 
points for each patch. Using a cubic spline passing through the 
key points, a desired grid-point distribution (clustering) is set 
up on the body surface. Then, by choosing an appropriate 
outer boundary, the grid for the flowfield calculation is 
generated by using an elliptic grid generator. More discussions 
can be found in Ref. 2. 

Results 

Results for the following five different configurations are 
presented to demonstrate the versatility and robustness of the 
code in handling a wide variety of nonlinear flows: 

Case I: Advanced tactical fighter configuration (Fig. 3). 

Case 2: Langley canard-wing fighter configuration (Fig. 8). 
Case 3: Isolated shuttle orbiter (Fig. 12). 

Case 4: Multibody configuration: shuttle orbiter with external 
tanks/solid rocket boosters (Fig. 17). 
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Fig. 4 ATF chord wise pressure distribution; Af w «1.6» a -1.24, 
Z=- 183 in. 



Fig. 5 ATF chordwise pressure distribution; \f„ =» 1.6, a = 1.24, 
2=245 in. 



The test conditions for each of these cases are summarized in 
Table 1. 

Case 1 

Figures 4 and 5 show the chordwise pressure distribution on 
the upper and lower surfaces at 60 Vo (183 in.) and 80 Vo (245 
in.) span stations, respectively. The results show that the pres- 
ent predictions are in very good agreement with Rockwell’s ex- 
perimental data. 

Figure 6 shows the comparison of overall forces and 
moments in terms of C L , C D , and C M . The full-potential 
results compare very well with NASA data. The drag calcula- 
tion C D includes the skin-friction and base drag. The com- 
puted viscous drag was for a fully turbulent condition at the 
test unit Reynolds number of 2 x 10 6 . The 48-deg wing sweep 
results of Fig. 6 correspond to a supersonic leading-edge con- 
dition of 3.3 deg. The Uft and drag coefficients from the pres- 
ent calculations for this fighter model, at different Mach 
numbers, are summarized in Table 2. The results are in ex- 
cellent agreement with experimental data. 

Figure 7 shows the grid and pressure contours for the same 
fighter geometry with a nacelle mounted on the undersurface 
of the wing. Only the exterior of the nacelle is modeled as part 
of the wing-body combination. At an axial marching station 
immediately preceding the inlet face, initial conditions are 
generated by interpolation from the flow field without the 
nacelle. The shock formed around the nacelle near the inlet 
face (see Fig. 7) is diffused at downstream stations. 

Case 2 

Figure 8 shows a fighter model tested at NASA Langley that 
has a canard and a fuselage-mounted flow-through nacelle. 
The actual computational geometry and the surface grid 
employed in this study are shown in Fig. 9. Computations 
were performed for this configuration at M m =2 and a = 4 
deg. Figure 10 shows results of cross-flow streamlines, surface 
pressures, pressure contours, and cross-flow velocity vectors 
at an axial station where the fuselage, wing, canard wake, and 
nacelle are all presented. The nodal singularity in pressure 
contour present at lower wing-body junction regions cor- 
responds to a saddle singularity of cross- flow streamlines, as 
shown in Fig. 10. Note the pressure match along the canard 
wake cut. The upper and lower center plane pressure contours 
at M m =2.0 and a = 4.0 deg are shown in Fig. 11. The bow 
shock, canopy shock, nacelle shock, and expansion wave are 
all nicely presented in this figure. 



Fig. 6 Comparison of measurement and full-potential prediction at 
1.6 for 48 deg sweep linear multipoint design wing. 


Table 2 Test cases for the advanced tactical fighter configuration 


— 

a 

5* 

4.5* 

5* 

4.5* 


M„ 

1.6* 

1.6 b 

I.4 b 

1.6 b 


A 

48* 

48* 

48* 

55* 

C L 

Code 

0.298 

0.3016 

0.3561 

0.29186 


Data 



0.342 

0.283 

Cp 

Code 

0.0482 

0.04916 

0.04117 

0.0404 


Data 

0.0457 

0.0493 

0.0426 

0.0396 


‘Without vertical tail. b With vertical tail. A = wing sweep angle. 


Case 3 

Figures 12-16 give the geometry, the gridding, and the cor- 
responding flowfield solutions of the isolated shuttle orbiter at 
- 1 .4, a = 0, and - 1 .94 deg. The chordwise pressures on 
the upper surface are shown in Fig. 13, and they compare very 
well with the experimental data. Figure 14 shows the cir- 
cumferential pressure distribution for the orbiter at x- 1200 
in.* It is noted that the pressure along the vertical tail and the 



Fig. 7 Pressure contour and grid of ATF with nacelle; M = 1.6, 
a = 5 deg at jc=375 In. 
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Fig. 8 Langley canard-wing fighter configuration. 



Fig. 9 Computational geometry and surface gridding for Langley 
fighter configuration. 



Fig. 10 Solution for Langley flgher configuration, Afo.^2.0, 
a = 4.0, x= 14.0; a) Streamline, b) Surface pressure coefficient, c) 
Pressure contour, d) Velocity vector. 


orbital maneuvering subsystem (OMS) pods are well 
predicted. 

Figure 15 shows the details of the orbiter geometry as 
modeled in this study. The OMS pod is clearly seen in Fig. 15. 
Figure 16 presents a series of isobar plots at different x loca- 
tions. The onset of the OMS pod shock formation is clearly 
seen. The OMS pod shock is formed around 1065 in., then 
grows, and finally hits the upper wing surface at approxi- 
mately x= 1090 in. The foot of the OMS pod shock moves fur- 
ther away from the fuselage for increasing x along the orbiter. 
The trace of the shock foot on the upper surface is also shown 



Fig. 11 Upper and lower centerplane pressure contour for Langley 
fighter configuration; Af* =■ 2.0, a = 4.0. 



Fig. 12 Shuttle orbiter configuration. 


in Fig. 16, and a comparison between the experimental shock 
and the numerical prediction is made. It should be mentioned 
that, since the present method is valid for supersonic flow with 
an embedded subsonic region, it allows one to treat the actual 
shuttle orbiter without having to make any modifications in 
geometry. 

Case 4 

Figure 17 schematically illustrates the multibody interaction 
problem of the shuttle orbiter in a mated configuration with 
the external tank (ET) and solid rocket boosters (SRBs) pres- 
ent. Figure 18 shows a perspective view of the complicated 
multibody problem as modeled by this full-potential code. It is 
found that the external tank has no effect on the upper wing 
surface and only a small effect on the lower wing surface of 
the shuttle orbiter. The high pressure present on the lower sur- 
face of the orbiter wing is caused by the aft attach struts that 
connect the orbiter to the external tank. The presence of the 
aft attach struts is modeled by a wedge blockage effect and the 
ET and SRBs are modeled by an elliptic cross-section exter- 
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Fig. 13 Shuttle orbiter upper surface pressure distribution. 



*(in| 


Fig. 14 Circumferential pressure distribution for the orbiter at 
x=* 1200 in., -1.4, « = 0. 




Fig. 16 Trace of OMS pod shock on the upper surface; M * => 1.4, 
a = 1.94. 


Fig. 20 Orbiter lower surface chordwisc pressure distribution; 
A/a, -1.4, a -0 deg. 
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nal tank. Figure 19 shows the gridding and- the isobar plot for 
this multibody problem at an axial station where the wedge 
blockage effect is present. The detached bow shock formed by 
the blockage is clearly shown in this figure. Chordwise 
pressure distribution of the orbiter lower surface at z = 0.34, 
with and without the blockage, is given in Fig. 20. The result 
with the blockage effect included shows a very good com- 
parison with the experimental data. 

Conclusion 

The main objective of this study is to illustrate the versatility 
and usefulness of a recently developed nonlinear aerodynamic 
prediction capability based on the full-potential equation. 
Results are shown for a variety of complex configurations, in- 
cluding a multibody problem. Comparison of results with 
available experimental data are in good agreement. The fully 
vectorized version of this code takes about 4 to 5 min of CPU 
time for analysis of typical fighter-like configurations on the 
CYBER 176 computer and about 25-30 s on the CRAY-XMP 
for a marching plane grid of 70 x 25. 
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A new finite-difference scheme has been developed to 
efficiently solve the Euler equations for three-dimensional 
inviscid supersonic flows with subsonic pockets. The tech- 
nique utilizes planar Gauss-Seidel relaxation in the march- 
ing direction and approximate factorization in the cross- 
flow plane. It is a unified formulation based on the unsteady 
Euler equations: an infinitely large* (‘infinitely small* re- 
ciprocal of) time step is used in parts of the flow-field 
where the component of velocity in the marching direction 
is supersonic — here the Gauss-Seidel sweeps are restricted 
to the forward direction only and the procedure reduces to 
simple space- marching; a finite time step is used in parts of 
the flow-field where the marching component of velocity is 
subsonic — here backward and forward Gauss-Seidel sweeps 
are employed to allow for upstream and downstream prop- 
agation of signals, and a time-asymptotic steady state is 
obtained. The discretization formulae are based on finite- 
volume implementations of high accuracy (up to third or- 
der) Total Variation Diminishing formulations. The fully 
general coordinate treatment used permits the use of ar- 
bitrary marching fronts (rather than just planes perpen- 
dicular to an axis, spherical fronts, etc.). Results are pre- 
sented for an analytically defined forbody, a twisted-cone 
Inlet spike, a realistic fighter configuration, and the Space 
Shuttle. 

1»Q In tro du ction 

For fully supersonic flows, an efficient strategy for ob- 
taining numerical solutions is to employ space-marching 
techniques. Realistic high speed flight vehicle configura- 
tions often give rise to subsonic pockets even though they fly 
at supersonic speeds. For such predominantly supersonic 
flows, a hybrid approach is suitable: a space marching tech- 
nique for the supersonic parts and a relaxation technique for 
the subsonic parts. Such a hybrid approach has been devel- 
oped for potential flows by Shankar, Szema et al 1,3,3 . For 
the Euler equations, however, the hybridization is conven- 
tionally achieved by coupling separate space-marching and 
time-marching codes, each with disparate grid systems, etc. 
Here, we present a unified approach for efficiently solving 
the Euler equations for three-dimensional supersonic flows 
with subsonic pockets. The aim is to develop an Euler 
solver as versatile as the potential flow solvers 1,3 ' 3 in treat- 
ing complex and realistic aircraft, space shuttle, and other 
types of flight vehicle configurations. By solving the Euler 
equations, however, we hope to be able to compute a wider 
range of flows with stronger shocks, rotational slip streams, 
etc. which void the irrotationality assumptions built into 
the potential flow simulations. 

* Member Technical Staff 
Computational Fluid Dynamics Department 
Member, AIAA 

Copyright © Aretricu lutlutc of Acro«a»tks iad 
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The new approach utilizes finite-volume implementa- 
tions of high accuracy (up to third order) Total Variation 
Diminishing (TVD) discretizations and are thus expected 
to be more accurate and reliable than other Euler space- 
marching and time-marching techniques based on central 
difference approximations. In contrast to these latter meth- 
ods, there are no parameters in our approach for fine-tuning 
numerical dissipation for every case. Numerical oscillations 
are, for the most part, eliminated by using TVD scheme 
based discretizations. 

The new approach is based on the unsteady Euler 
equations. However, in the supersonic parts of the flow 
(where the velocity component normal to the cross-flow 
plane that identifies the local marching direction is super- 
sonic) an ‘infinitely large* time step (which implies an ‘in- 
finitely small’ reciprocal of time step) is employed. This 
makes the transient terms of the discretized unsteady equa- 
tions vanish. In subsonic parts of the flow, a finite time step 
is employed and the steady-state is approached as a time- 
asymptote. 

The new solution approach is based on a planar Gauss- 
Seidel relaxation method coupled to approximate factor- 
ization in the cross-flow plane. In supersonic parts of the 
flow-field, the Gauss-Seidel method is restricted to forward 
sweeps and thus the solution procedure reduces to a sim- 
ple marching technique. In subsonic parts, both forward 
and backward sweeps are used along with the finite time 
steps mentioned earlier. Stability of such an approach is 
guaranteed by the diagonal dominance resulting from us- 
ing TVD discretizations in the marching direction in the 
transonic parts of the flow-field. This is a crucial difference 
between conventional hybrid Euler solvers and the new ap- 
proach. In conventional approaches, space-marching and 
time-marching techniques must be applied in overlapping 
regions for stability. In the new unified approach, there is 
no need for overlap. 

In the following sections, we describe the new method 
in detail. We first cast the equations in finite-volume dis- 
crete conservation law form. Then we explain how the vol- 
ume and the metrics are evaluated. This essentially com- 
pletes the treatment of geometry and we proceed next to 
the details of the algorithm. TVD discretizations are ex- 
plained first. Then the marching/ relaxation procedure is 
described. This covers the use of approximate factorization 
in the cross-flow plane, the reduction of the Gauss-Seidel 
procedure to a marching procedure in supersonic zones, etc. 
The boundary point treatment is also explained briefly. 

In the results section, calculations for an analytically 
defined forebody are presented first to illustrate some fea- 
tures of the new algorithm. Results for many conical flow 
cases have also been obtained but are presented elsewhere 4 . 
Next, results for a twisted-cone inlet spike are shown. Then, 



results are presented for a realistic fighter aircraft configu- 
ration with fuselage, canopy, wing, nacelle and vertical tail. 
Finally, we conclude by presenting results for the Space 
Shuttle Orbiter configuration. 


2.Q.. Finite.vplumc Framework 


In this section, we describe the finite-volume frame- 
work chosen to implement the algorithm. We start by in- 
troducing the semi-discrete conservation law form and asso- 
ciating it with a finite-volume formulation of the geometry. 
Then we provide detailed formulae for the evaluation of the 
cell volume and cell-face normals. 

21 S« mi‘digc retc Cpn ggyatfon Law 


We begin with the conservation law form of the un- 
steady Euler equations in the Cartesian coordinates x, y, r, 
and time t 


Qt + -S* + + G m = 0 


(2.1a) 


where the dependent variable vector <J, and the fluxes E , 
F % and G are given by 


F = 



/ « > 


/(e + p)u\ 
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pu 

Q = 

pu 


pu 2 + p 


pv 


pvu 


\pu)j 


V pwu / 

/(« + pM 


/(e + p)uO 

[ p» 


pw 

puv 

,G = 

pUW 

p* 

+p 


pvw 

V pwv / 


l piv 2 + p J 


( 2 . 1 b) 


In the above, pressure is p, density is p t Cartesian x, y, z 
velocity components are u,v, tn, and the total energy per 
unit volume is t computed from e = p/( 7 — 1) + p(u 7 + v 1 + 
u> 3 )/2. 

Assuming a time invariant grid, under the transforma- 
tion of coordinates implied by 


€ = {(*. If.*). n = ti{x,y,z), 


( 2 . 2 ) 


Eq. 2.1 can be recast into the conservation form given by 


where 


Q r + + ?n + Gf — 0 


»-! . 

E=^E+^fF+^G 

J J J 
F = -J-E+ ~j~E+ ^—G 

G=$jE+&F+^G 


(2.3a) 


(2.36) 


where, in turn, J is the Jacobian of the transformation 

J -*(*.*, *)/*(*.*.*) (2.3c) 

Associating the subscripts ;,*,/ with the direc- 

tions, a numerical approximation to Eq. 2.3a may be ex- 
pressed in the semi-discrete conservation law form given by 

(Qj,*,/)r + (-®j+t/3,JM “ £>-l/3,fc,l) 

+ (^y,*+ 1/3,1 - */,*- 1/3,0 (2.4) 

A .A. 

+ (^j.fc.l+1/3 ” Gj, *,1-1/3) = 0 

where P t P, G are numerical or representative fluxes at the 
bounding sides of the cell for which discrete conservation is 
considered, and Qy.fc.i is the representative conserved quan- 
tity (the numerical approximation to Q) considered con- 
veniently to be the centroidal or cell-average value. The 
half-integer subscripts denote cell sides and the integer sub- 
scripts the cell itself or its centroid. In Fig. 1, the eight 
vertices of one computational hexahedral cell are identified 
by numerals 1 through 8. These must be associated with 
the appropriate /, *, l triplets: 

1=;- 1/2,*— 1/2,/- 1/2 
2 = y + 1/2, * — 1/2, / — 172 

3hj- 1/2,* + 1/2,/- 1/2 

4 = 1/2,*- 1/2,/+ 1/2 

5 = y+ 1/2,*+ 1/2,/- 1/2 1 ■ } 

6 =i- 1/2,* + 1/2,/+ 1/2 

7 = y + 1/2,* — 1/2,/+ 1/2 

8 = y + 1/2,* + 1/2,/+ 1/2 

In the following, subscripts easily understood by implica- 
tion will be dropped for brevity. 

The semi-discrete conservation law given by Eq. 2.4 
may be regarded as representing a finite-volume discretiza- 
tion if the following associations are made: 



Fig. 1 Computational finite-volume cell 
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Qi.kj = QV, ik ., (2.6a) 

where V is the volume of the cell under consideration; 


(kiA = 

\ J ' i± i/3 

- 1/2, l - 1/2), (* + 1/2, 1-1/2), 

(* + 1/2, / + 1/2), (* - 1/2, / + l/2)} y±l/a , 

/ Vx,y,x \ 

V J /fc± 1/2 

n *,y,*{0 - 1/2, 1-1/2), (; -1/2,1+ 1/2), 

U + 1/2, 1 + 1/2), (j + 1/2, 1 - 1/2 )}*±i/2 , 

/ = 

V J /J±l/2 

- 1/2, * - 1/2), (; + 1/2, k - 1/2), , 

(j + 1/2, k + 1/2), (; - 1/2, fc + l/2)}, ±1/3 . 

( 2 . 66 ) 


In the above, n Xt ,., are the x, y,z components of the repre- 
sentative normals to the surface formed by the four points 
a, 6, c, d implied in n X( , iX (a, 6, c, d). Four points do not nec- 
essarily lie in one plane and therefore the components n Xt y iS 
refer to representative values for an equivalent single plane. 

The evaluation of the volume and metrics (cell-face 
normals) are now presented in the following subsections. 
The evaluation of the representative flux is presented in 
the next major section. 


2.2 Computation of Cell Volume 


First, the volume of a tetrahedron denoted by its ver- 
tices a, 6, c, d is evaluated from 


V“ t {a,b,c,d) = 

kaMZe - Zi) ~ Vc{Zb ~ *d) + Vd{Zb ~ 2 C )1 
-2k[y 0 (2c - zd) - y c (z a - id) + y<f(z« - z c )l (2.7) 

+2c[y a (2k - Zd) - y*(Za - Zd) + Vd[z a “ «t.)] 

-Zd\y a {zb - Z e ) - y h (z a - z c ) + y c [z a - Z6)]|/6.0 


Then, referring to Fig. 1 again, the volume of the hexahe- 
dron is computed as a sum of the six individual tetrahedra 
that constitute it. 

V = V tet [\,2,S,7) + V tet {l,7 f S,S) 

+ V* e *(l,3,5,8) + K‘*‘(l, 3,8,6) (2.8) 

+ K'** (1,7, 8,6) + V* e *(l, 4,7,6) 


It is of interest to note that such a formula will result in the 
proper evaluation of volume even when some of the faces of 
the hexahedron collapse to a line or a point. 

2.3 Computation of Cell- Face Metrics (Normals) 

In Eq. 2.6b, cell-face normals were introduced. Each 
cell face is identified by four vertices not all of which are 
necessarily on a single plane (three points being sufficient 
for defining a plane). In our approach, we allow for this 
and also for some of the faces to collapse to an edge or even 
a point. Computationally, we will always identify a face 
by its four vertices a, 6, c, d expressed in the j t k , l subscript 


system. Physically, some or all of the four vertices may lie 
at the same x, y, z location. 

The cell-face normals are evaluated as 

n x (a, 6,c,d) = [dy^dz^ - dy^dz ^)/ 2 
+ (dy<udz a d — dyaddzdc)/! 
n f (a, 6, c, d) = (dz^dx e b - dz^dx 6«)/2 
+ ( dZdcdXad dZaddXdc)/! 
n,(a, 6, c, d) = (dri*dy c 6 - dx^dy^)/! 

+ (dx<udy a d - dx ad dydc)l 2 

where 

dSj 3 “ Si — S a 

where, in turn, s corresponds to x or y or z, and 1 and 2 
correspond to a or 6 or c or d. The first term in each of 
the definitions is respectively one half of the x, y, or z com- 
ponent of the cross product of the vector from a to 6 with 
the vector from 6 to c. The second term in each definition 
is correspondingly one half of the x, y, or z component of 
the cross product of the vector from c to d with the vector 
from d to a. A cross product of two vectors lies along the 
direction of the normal to the two vectors. In the present 
situation, the two vectors are connected. Therefore, such a 
normal also defines the direction of the normal to the plane 
containing the two vectors. One half of the cross product of 
connected vectors also has a magnitude equal to the area 
of the three dimensional triangular planar shape defined 
by the two vectors. Thus, while {*, y ,*/d t f; x , yil /J, U.y.x/J 
define the x, y, x components of the normals (not unit nor- 
mals) to the local tangent plane to the constant sur- 
faces respectively, the associated quantities (n x?yiX )y,*,i de- 
fine the components of the normals to the local constant 
;, k t l planes. 


(2.9a) 


(2.96) 


3.0 TVD Discretization 

In the last section, the numerical or representative 
fluxes E,F t G were introduced. These flux_es_are so named 
because they approximate the real fluxes E } F, G to the re- 
quired order of accuracy. The actual fluxes appearing in 
the governing partial differential equations depend on the 
metrics and correspondingly, we 

allow the numerical fluxes to depend on the numerical met- 
rics (the cell-face normals). In the last section we did point 
out the link between the metrics and the components of the 

cell-normals, but the numerical flux was not defined there. 
The latter task is the subject of this section. 

We employ an upwind-biased scheme in our approach 
in such a fashion as to essentially eliminate numerical or 
spurious (unphysical) oscillations while, at the same time, 
achieving high accuracy. In order to describe this type 
of discretization, we first mention the underlying upwind 
scheme used in terms of the corresponding approximate 
Riemann Solver, and then expand upon how high accu- 
racy and the TVD property’ are built in. More details on 
these and related topics may be found in Refs. 5 and 6 and 
in references cited therein. 


3 



3.1 Roe’s Approximate Riemann Solver 

The Riemann Solver is a mechanism to divide the flux 
difference between neighboring states (between Q m and 
Qm+i for e.g.) into component parts associated with each 
wave field. These can in turn be divided into those that 
correspond to positive and negative wave speeds. When we 
compute the numerical flux at the cell face at m+ 1/2 using 
various combinations of fluxes and positive and negative 
flux differences, in the finite-volume formulation, we will 
only use the cell-face normals defined at m + 1/2 in all the 
terms contributing to that representative flux. The actual 
fluxes E , F , G, when evaluated with the metrics equated to 
cell-face normals, can all be written in the same functional 
form given by 


At each cell face, the positive and negative projections 
of the eigenvalues may be defined by 


„■+ (^rn+I/2 ^ l^m+l/aD 

A “ 2 


,»* = l.-,5 


(3.5) 


In order to help Roe’s Riemann Solver avoid expansion 
shocks, only at sonic rarefactions (A'(<? m , W m +i/ 2 ) < 0 < 
A’(Qm+i, JVm+i/a)), the corresponding positive and nega- 
tive projections are redefined as 


V± — l 1 * 

.rn+l/2 — A m+ 1/2 

± (A’fQm+i, JVm+1/2) ~ A t (Qmt^frH-l/2)) 


= /»,»«. n„n.) = f(Q t N) (3.1) 

where the appropriate values of n,,nj,,n x are used and N 
denotes the set of those normals. Using such notation, it is 
possible to present the necessary algebra very concisely. 

Let us first denote the Jacobian matrix of the flux / 
with respect to the dependent variables Q by df jdQ. This 
Jacobian can also be called the coefficient matrix. Let us 
denote the eigenvalues of the coefficient matrix by A* and 
the corresponding left and right eigenvectors by C and r*, 
respectively. The matrix formed by the left eigenvectors. as 
its rows is then called the left eigenvector matrix L and the 
matrix of right eigenvectors comprising right eigenvectors 
as its columns is R . For our purposes, we choose an or- 
thonormal set of left and right eigenvectors which implies 
that LR — RL = /, the identity matrix. In the above, the 
superscript t has been used to denote the association of the 
«*th eigenvalue with its corresponding eigenvector. Each 
eigenvalue is also associated with its own wave field. 

The underlying upwind scheme is based upon Roe’s ap- 
proximate Riemann solver 7 . In this approach, cell interface 
values of density, velocities, and enthalpy (A = 7 p/((i - 
l)p) + (u 3 + t/ 3 4- v> 2 )/ 2) are computed using a special av- 
eraging procedure: 


Pm+1/2 


(tl J V l w)m+,/ 2 = 


__ Pm+ 1 \J Pm + 1 + Pm\/pm 
y/Pm+l + \fPrn 
(u»V)«Om+ly/Pfn+l 


n+1/2 = 


y/Pm+i + y/p^ 
_ ^m+1 >/Pm+\ + h my /p m 
y/Pm+l + yfp^x 


(3.2) 

where m = ) or k or /. From the above, the speed of sound 
can be computed from 

Cm+i/a = \J {Am+ 1 /a - («* + U 3 + u> 2 )/2}(7 - 1) (3.3) 


Knowing (u, v, w, e) m+ ,/ 3 , the eigenvalues and orthonormal 
set of left and right eigenvectors corresponding to a cell face 
can be computed. These may be denoted by 

^m+l/2 = ^m+l/2(^?"»+l/2» ^m+l/2)» 

^m+1/2 = Cn+l/2Wm+l/2i^m+l/2)> (3.4) 

r m+l/2 = r m+I/2(Cm+l/2i Af m+1 / 3 ). 


For the sake of completeness, detailed formulae for the 
eigenvalues and the eigenvector matrices are now presented. 
Defining the contravariant velocity by 

17= n x u + n y v-f n x w , (3.7) 

the eigenvalues are given by 

A 1 = U - cyjn 2 x + n 3 + n 3 
A 3,3,4 = U (3.8) 

A^tZ + c^/nl + nJ + nJ 

Defining 

n x , v .x = ni, y ,«/^/n 3 + n 3 + n 3 (3.9) 

and 

t? = (u 3 + t/ 3 + u> 3 )/2 , (3.10) 

the left and the right eigenvector matrices are given in Ta- 
ble 3.1 and Table 3.2 respectively. 

3.2 High- Accuracy TVD Schemes 

We can construct upwind-biased schemes of varying 
accuracies using the basic ingredients given in the last sub- 
section. Here, we present a family of schemes based on 
the preprocessing approach 4 . Let us now define some con- 
venient variables as an intermediate step before defining 
the numerical flux corresponding to a high-accuracy TVD 
scheme. First, we define a parameters which are a measure 
of the change in dependent variables across the correspond- 
ing wave family and therefore are a measure of the slope 
between neighboring states. In the following, the super- 
script i corresponds, as usual, to the i-th eigenvalue and 
i-th eigenvector. The subscripts 1, 2, and 3 are just la- 
bels to differentiate between the three different types of a 
parameters. 

a l,m+l/2 = ^m+l/2(*?m “ Qm-l)i 

Q 2,m+l/2 = 4n+l/2(Qm+l " Qm ) i (3.11) 

a 3,m+l/2 = Ci+l/2(*?™ + 2 _ Qm+ 1)* 
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Table 3.1 The left eigenvector matrix L 


Next, we define the slope-limited values given by 

^l.m+l/3 = m ^ nm0 ^l a l,m+l/3»^ a 3,m+l/3]» 

a 3.m+l/3 = n ^ nm0 d[ o 3,m+l/3»^0 f l,m+l/3]* . . 

. . t 312 ' 

a 3.m+l/3 — n ^ nmoc ^[ Q 3,m+l/3» ^ a 3,m+l/3]» 

“3,m+l/3 = 

In the above, the compression parameter 6 is to be taken as 
the following function of the accuracy parameter <j> which 
is explained shortly. 

6 =r^ < 313 > 

The “minmod* slope-limiter operator is 
minmod[x,y] = sign(x) max[0,min{|x|,ysign(x)}] (3.14) 

In Eq. 2.4, we introduced numerical fluxes & % P,G, 
Basedjon the concise notation of using / to represent ei- 
ther E or F or G^ let us use / to denote the numerical 
fluxes S or P or G. We can then write down a family of 
TVD schemes as follows in terms of the previously defined 
a parameters (with the subscript m + 1/2 dropped from 
these for convenience): 


^m+1/3 — j [/W m + l» W m+ i^a) + /(Qmt ^m + 1/2)] 
~ 2 S (^m+l/3 “ \n+l/3) a 3 r m+l/3 


— f(Qmt ^rn+l/3) + ^m+ 1 ^ 3 a 3 r m+l /3 

t 

= /(Qm+1 j ^m + I/3) ~ ^m+I/3 a 2 r m+l/3 

(3.16) 


The remaining terms on the right hand side of Eq. 3.15 
define correction terms that upgrade the accuracy. For use 
in the next subsection, we define 

^/m+l/3 = ^m+l/3 a 3,m+l/3 f m+l/3 * (3-17) 


It is interesting to note that in all the above formulae used 
to define the numerical flux at m + 1/2, the eigenvectors 
and eigenvalues are only necessary at the corresponding 
cell interface. Therefore, the only geometry information 
used corresponds to the cell-face normals at m + 1/2. The 
solution variables Q are sampled between the centroidal 
points m— l,m,m+l,m+2 when the various a parameters 
are defined. 


/m+t /3 = ^m+ 1/2 

,S ^( i + igi , 1 ~ ^ g* \ > i+ t 

+ 2^1 4 "*■ 4 a * Mm+i/a r m+i/a 

»■ ' ' 

(3.15) 

The first term on the right hand side of Eq. 3.15 defines 
a first-order numerical flux and is constructed from 


The parameter <f> defines schemes of varying accuracy. 

The notations at* and 5 have been used to define slope- 
limited values of the a parameters. If we replace these by 
their unlimited values, we obtain schemes whose truncation 
error in one-dimensional steady-state problems on uniform 
grids is given by 

• (3,.) 


5 




\* e 

- + - n x u 

(c 7-1 

- h y v - h x w fs/l 

+n x ti — n x w 

7"* 

-n„u + h x v 

-h z v + h v w 

t? c 

- + 7 + n x u 

c 7-1 

+ fiyV + n,ty /\/2 

[i]/v /2 

c 

e 

c 

- 1 />/* 
c 


u - , - 
— riy 4 - 

«. . 
-n,-n v 

u A 

7”* 


— + n x j /y/2 


v. 

c n » 

^ n, + h x 

^n x -h. 


J +n y ]/v / 2 

(i -*■]/* 

w . 

- n x 

w . 

-n, 

c 

w . 

-rij + fly 

c y 




Table 3.2 The right eigenvector matrix R 


Here, the truncation error refers to the difference between 
the centroidal value of the numerical solution and the av- 
erage value of the exact solution in that cell. The choice 
of (f> = 1/3 results in a TVD scheme based on an under- 
lying third-order scheme. The choice of <j> = -1 results 
in a TVD scheme based on the fully upwind second-order 
accurate formulation. FVomm’s scheme arises when $ = 0 . 

3.3 TVD Schemes and Diagonal Dominance 

In the next section, a procedure is presented to solve 
the finite difference equations resulting from the TVD dis- 
cretization of the space differencing terms. In supersonic 
zones, the method reduces to a simple marching scheme, 
while in subsonic zones it becomes a relaxation approach 
and both forward and backward sweeps are employed along 
the marching direction. In order for such a relaxation ap- 
proach to be stable, a sufficient condition is the diagonal 
dominance of the underlying finite difference scheme. This 
diagonal dominance can be shown to exist for TVD dis- 
cretizations. For more details, the reader is referred to 
Ref. 8 . 


4.0 The S olution Procedure 


We begin this section by considering an implicit time 
discretization coupled with the TVD space discretization 
discussed earlier in terms of the corresponding numerical 
flux terms. 


Q n+I - Q n 
At 


+ 


+ 


+ 


(^y+i/a-^y-i/a) 
(ft+i/a -ft-i/a) 

^ \ n+l 

Gi+i/a =0 


(4.1) 


Here, n is the index in time and At is the time step. In what 
follows, we will consider the linearization of the above non- 
linear set of finite difference equations. Then we will sim- 
plify the algebraic solution procedure by approximately fac- 
torizing the implicit operator in the cross-flow plane (which 
is a plane in only computational coordinates — constant j 
plane). The marching direction is along We will further 


specialize the scheme for the two cases of supersonic and 
subsonic velocity components in the marching direction. 

4.1 L inear is at io n 

Let us linearize Eq. 4.1 about a known state Q = q 9 
using a Newton procedure to obtain a better approximation 
q‘+ l to Q n+1 . Here, s is a subiteration index. Defining 

A 'q = q‘ +l -q‘ 

Ay£ = Ej+if 7 - Ej-lf2 
A kF = ^Jt+1/2 - -Rfe-l/2 

AiG = G/+1/2 - £?f_i/2 i 
we can describe the Newton procedure by 


£ : /+^(a,£ + a*/ , + a,g) 


a *9 = 


— - (?") + Ay £(?') + A k F(q‘) + A,G{q') 


( 4 . 3 ) 


We next simplify the left hand side by defining an ap- 
proximation to the true linearization. Towards this goal, 
we consider only a first-order accurate scheme (based on 
the first-order numerical flux A) for the left hand side while 
we include the full high-accuracy scheme on the right hand 
side. Even so, when the subiterations converge, the right 
hand side is satisfied to the desired degree. Next we assume 
that the eigenvalues and eigenvectors are not functions of 
q. Finally, we observe that 


y! ^m+l/2 “ Jbn— 1/3 = 


E 1/2 + E E^m+ 1/2 


(4.4) 


because, in expanding Eq. 4.4 using Eq. 3.16, we find that 

y [(n*)m+l/3 (tt*)m-l/2l = ® 

E K n v)m+l/2 ” ( n v)m — l/a] = 0 (4 5) 

m=y,fc,Z 
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yi [( n *)m+l/3 ” (n#)m— l/i] _ 0 
m=},k,l 

when the cell-face normals are evaluated using the formulae 
given in Eq. 2.9. Using the above, Eq. 4.3 is simplified to 



+ B* +l/3 (A*7 fc+1 - A'**) 

+C+ l/3 (A'* - A'*- 1) + C,- l/a (A*g /+l - A'*) 


= Right Hand Side of Eq. 4.3 


where 


(4.6) 

A f± 1/3 = • R /±J/3 A *± l /2 i j±»/a 

B k± 1/3 = £*±l/3 A *±l/ 3 £ *±l/2 

(4.7) 



Here, 

A ± = (A + |A|)/2 

(4.8) 


in which A is the diagonal matrix whose diagonal elements 
are A* and |A| is the diagonal matrix whose diagonal ele- 
ments are |A'|. 


4.2 Planar Gauss-Seidel Relaxation 

Even after the many simplifications leading to Eq. 4.6, 
it is obvious that more algebraic simplification is needed be- 
fore a computationally feasible and efficient solution proce- 
dure is obtained. This is because Eq. 4.6 signifies a system 
of equations which links every point /, k t l with its six neigh- 
bors y+l,/-l f *+l f A-l,/+l l /-lin such a fashion that 
the left hand side of Eq. 4.6, when considered for all grid 
points, is a huge (even though sparse) matrix whose band- 
width is also very large. Of course, for supersonic flows, 
a fully upwind difference approximation arises in the j di- 
rection and the dimensionality is reduced because the left 
hand side does not link j with j 4* 1. However, with our 
expressed aim of developing a method for subsonic pockets 
also, it is necessary to consider the case when j is linked 
with both its neighbors j — 1 and j 4-1. In such a case, a di- 
rect Gaussian elimination procedure for the matrix system 
of equations would be unacceptably expensive. Therefore, 
instead of a direct elimination procedure, we seek to ob- 
tain an efficient relaxation solution to Eq. 4.6. We choose 
a planar Gauss-Seidel procedure by retaining all terms of 
the left hand side except the off-diagonal terms in j (those 
terms that multiply A*$y±i). That such a procedure will 
be stable for TVD discretisations was discussed in Section 

3.3 and in the references cited therein. 

The planar Gauss-Seidel procedure can be written as 
+ 7 B *+i/a( A *9t+i “ 

v * t a n\ 


+ 7<£ l/3 (A'? ( -A'*-0 

+ 7<T +l/3 (A* 9l+1 -A* 9l ) 

= y [Right Hand Side of Eq. 4.3] . 

Denoting 

i=^4-yAt_ 1/a 4- yj47 +1 / a , (4.10) 

we can rewrite Eq. 4.9 as 

[/+ {s* + _ 1/3 a*_ 1/3 + Br +l/s A* +l/J 

+c j-i/a A <-»/a + c i+i/a A <+i/a}] A * 9 

= —A” 1 [Right Hand Side of Eq. 4.3] . 

(4.11) 

Of course, when the relaxation cycles denoted by super- 
script s converge to the desired extent, A *q = 0, and the 
full accurate formulae of the right hand side will be satisfied 
to a corresponding degree. 

4.3 Approximate Factorisation in the Plane 

While Eq. 4.11 defines an algebraic set of equations 
whose dimensionality is one order less than that of Eq. 4.9, 
it is still too huge to be tackled by an elimination algorithm. 
Therefore, we will now further reduce the dimensionality by 
approximately factorizing the left hand side of Eq. 4.11 to 
result in 

I + A* +1 /j j] 

/ + 7 ^ |c,t t/ jA)-i /2 + C ( " l/3 A| + i/j}] A'g 
= [Right Hand Side of Eq. 4.3] . 

(4.12) 

The actual sequence of steps to solve Eq. 4.12 can be chosen 
so that j4 “ 1 need not actually computed and only A is 
needed. For this purpose, we solve, in order, the equations 

[^+7 + B *+1/2 A *+>/j}] 5 (4 13 ) 

= —[Right Hand Side of Eq. 4.3] 

and 

[*' + V + C l+»/a A <+»/a}] (4.136) 

= Aq . 

with q being a temporary storage variable. 

Let us summarize the solution procedure developed in 
Eq. 4.13 for just one constant j plane. Equation 4.13a 
must be solved for all A:*varying lines (for all /). Then 
Eq. 4.13b must be solved for all /-varying lines (for all values 
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of k ). However, each fc- varying or /-varying line is associ- 
ated with only a one-dimensional block-tridiagonal system 
of algebraic equations whose block matrices are 5 x 5. These 
two steps only constitute one cycle of the Gauss-Seidel it- 
erations and that too for only one constant j plane. The 
planar Gauss-Seidel procedure requires that one constant j 
plane is updated at a time. When the neighboring j plane 
is updated next, the latest available values of the update 
variables q are used in the right and left hand sides. The j 
sweep strategy will be specialised for supersonic and sub- 
sonic how regions in what follows. 

4.4 Programming Notes 

We store grid information at two planes (grid-planes 2 
and 3) which describe the j boundaries (;' - 1/2, j -f 1/2) 
of one plane of cells. Let the centroids of these cells be 
denoted as solution-plane 3. Array storage is provided for 
dependent variable planes (solution-planes) 1,2, 3, 4,5. As 
the solution is marched, the contents of the grid-plane and 
solution-plane arrays are updated by replacing them with 
neighboring values or by the planar Gauss-Seidel algorithm. 

Let us consider the very first marching sweep now. We 
begin by initializing the two grid planes and the dependent 
variables at solution-planes 1 and 2. We are interested in 
updating the solution at plane 3. We first set the solution 
at planes 3, 4 and 5 to be equal to the values at solution- 
plane 3. After one or more subiterations for solution-plane 
3, we shift our attention to the next j -plane. Grid-plane 2 
is replaced with the contents of grid-plane 3. Grid-plane 3 
nodal values are stored on auxiliary storage for later use. 
New values for grid* plane 3 are generated by grid genera- 
tion procedures or read in from auxiliary storage initialized 
previously. Similarly, solution-plane 3 is saved on auxiliary 
storage for subsequent processing. Solution-plane 1 is re- 
placed by contents of solution-plane 2, and plane 2 is then 
replaced by contents of plane 3. Solution-planes 4 and 5 
are set to the values at plane 3 and the marching proceeds. 

If more than one subiteration is to be performed in 
the first marching sweep, the grid information is not up- 
dated for the subsequent subiterations. Solution-planes 4 
and 5 are reset to values at solution-plane 3 after the pre- 
vious subiteration and the next subiteration is processed. 
Solution-plane 3 values are not set to solution-plane 2 val- 
ues for the second and subsequent subiterations. 

For fully supersonic flows, a fully-upwind, not-flux- 
limited differencing scheme is used. Thus, the values set for 
solution-planes 4 and 5 are actually not used at all. For- 
ward marching is enough. Even first-order upwind scheme 
in the j - -direction and one subiteration per marching plane 
are also often enough. A small value is input for the recipro- 
cal of time step. Accuracy of approximate factorization for 
any time step size is maintained due to reasonable marching 
step size (distance between j grid planes). 

Subsonic regions could develop as a result of gradual 
compression (for e.g., around canopies) or abrupt transition 
through a shock wave (for e.g., in front of a blunt nosed ob- 
ject in an oncoming supersonic flow). In such regions, a 
larger value is chosen for the reciprocal of time step. The 
solution is marched forward using one or more (usually a 


maximum of two) subiterations by conforming to the pro- 
cedure outlined above for the first marching sweep. Then, a 
backward marching sweep (or even another forward march- 
ing sweep) is performed. For all sweeps (forward or back- 
ward) after the first, solution planes are filled with previous 
sweep solution values before updating using subiterations. 
Shifted replacements of solution-plane values of dependent 
variables are not used. For subsonic regions (subsonic pock- 
ets in supersonic flow), a TVD formulation of the desired 
accuracy is used enabling even strong shocks to be captured 
routinely. 

For very small pockets of subsonic flow caused by grad- 
ual compression, one forward sweep followed by one reverse 
sweep is enough. Even the reverse sweep is usually redun- 
dant in this case. For larger subsonic zones, a few (tens) of 
sweeps usually suffice. Residues are monitored for conver- 
gence. 

5.0 Boundary Point Treatment 

Only an outline of the boundary point treatment will 
be presented here due to lack of space. The boundary 
method used is fully compatible with the interior point dif- 
ferencing. It is based on considering a Riemann Initial and 
Boundary Value Problem at the boundary to construct the 
boundary point discretization. In this, it is similar in spirit 
to the correspondence between interior point discretization 
and the Riemann Initial Value Problem. The implemen- 
tation is specifically tailored to approximately factored im- 
plicit schemes. Linear boundary conditions (such as surface 
tangency) are exactly satisfied after every marching step. 
Corner points are also properly treated. More details on the 
new treatment used here including theoretical background 
and implementation details for explicit and other implicit 
methods are available in Ref. 9. A brief description of the 
type of boundary condition techniques used here can also be 
found in Ref. 5. Another approach to boundary condition 
procedures which can be applied to implicit schemes for the 
Euler equations is presented in Ref. 10. The importance 
of proper and accurate boundary condition procedures is 
demonstrated in Ref. 4. 

6.0 COMPUTATIONAL EXAMPLES 

The preceding sections have described an Euler March- 
ing Technique for Accurate Computations (EMTAC) and 
we now present many computational results obtained using 
the EMTAC code. The first set of results are for an an- 
alytically defined forebody geometry and these results are 
compared with experimental data. The next case consid- 
ered is the supersonic flow over a twisted-cone spike of a hy- 
pothetical aircraft inlet and the results are compared with 
numerical results obtained using a full potential marching 
code. The third set of results are for a realistic fighter con- 
figuration and once again most of the comparisons for this 
case are with the full potential marching code. The last set 
of results are for a Shuttle Orbiter configuration and the 
numerical results are compared with experimental data for 
this case. 
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6.1 Analytic Forebody 

Figure 2a shows the developed cross-section of a fore- 
body geometry reported in Ref. 11. The surface pressure 
distributions in the axial direction on the upper (9 = 0°, 
leeward side) and lower (0 = 180°, windward side) planes 
of symmetry at = 2.5, a = 0° are given in Fig. 2b. 
The grid and circumferential pressure distribution on the 
body surface at x/t = 0.22 and x/t - 0.34 for the same 
free-stream conditions are presented in Figs. 2c and 2d re- 
spectively. Figure 2e shows the circumferential pressure 



2.235 J 
(0.880) 




Fig. 2b Axial pressure distribution 


distribution on the same geometry for — 1.7, a — -5° 
at x/t = 0.278. It is noted that a small subsonic pocket 
develops, for this second case, on the lee side and two 
global marching sweeps are enough for the present numeri- 
cal method to give a very good converged solution. The ex- 
perimental data 11 are also presented in these figures. The 
comparisons show that the present numerical predictions 
are in excellent agreement with experimental data. 



Fig. 2c Circumferential pressure distribution 
at x/t = 0.22 



THETA 


Fig. 2d Circumferential pressure distribution 
at'x/£ = 0.34 










Fig. 2e Circumferential pressure distribution 
at Mqo = 1.7, x/^ = 0.34 


6.2 Twisted-Cone Inlet Spike 

Figure 3a presents the geometry of a twisted-cone in- 
let spike. At Moo = 2.5, the pressure contours at various x 
locations are given in Fig. 3b. The circumferential pressure 
at x ss 40 is compared with results obtained using a full 
potential solver (SIMP, described in Refs. 1-3) in Fig. 3c. 
As expected, the full potential method predicts a higher 
pressure on the upper surface where a strong nonisentropic 
shock is formed for the case considered. This shows the 
importance of using an Euler solver rather than a full po- 
tential solver for supersonic flow computations which must 
capture strong shock waves. 


6.3 Realistic Fighter Configuration 

Figure 4a shows the geometry and surface gridding of a 
realistic fighter-type configuration which includes a nacelle 
and a vertical tail. To illustrate the important features of 
the present analysis method, results have been obtained 
for the free-stream condition Mpo = 1.6, a = 4.94°. The 
results are compared with those obtained using the SIMP 
full potential solver. Figures 4b and 4c present the surface 
pressure at the upper and lower symmetry plane. The re- 
sults show the excellent agreement between the predictions 
of these two codes. Circumferential pressure distributions 
and pressure contours at two different x locations which in- 
clude the nacell, vertical-tail, wake and wing are presented 
in Figs. 4d-e. The comparison shows very good agreement 
except at the lower surface of the wing in the vicinity of the 
wake region. A higher pressure is predicted by the Euler 
(EMTAC) code. It is also noted that the wake treatment 
in both methods provides the correct zero pressure jump 
across the wake. 

Table 6.1 shows the comparison of overall forces in 
terms of Cl,Cd and Ci/Cp. The drag calculation in- 
cludes skin friction drag estimated using a boundary layer 
technique and an estimate of the base drag. Both the full 
potential and Euler results agree very well with Rockwell 
experimental data with the Euler results being closer to the 
data. 
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SIMP 

EMTAC 

DATA 

Cl 

0.30588 

0.3017 

0.303 

Cd 

0.032458 + 0.013 

0.03433 + 0.013 

0.0475 


= 0.045458 

= 0.04733 


Cl/Cd 

6.72 

6.38 

6.42 


Table 6.1 Comparison of Potential, Euler and 
and experimental data for Cl % CdiCl/Cd 


6.4 Space Shuttle Orbiter 



Figures 5a*5g give the geometry, gridding and corre- 
sponding flow-field solutions for an isolated Space Shuttle 
Orbiter flying at Afo© = 1.4, a = 0°. The EMTAC code 
is applied to compute the flow field about the entire or- 
biter, from nose to tail. Multiple (uni- or bi-directional 
sweeps are used in the nose region to capture the detached 
bow shock and the subsonic region behind it. After this 
subsonic region transitions by expansion, over the shoul- 
der region of the nose, into a supersonic flow-field, a simple 
forward-marching technique is employed. Multiple relax- 
ation sweeps are also used in the canopy and OHMS pod 
regions to compute the locally subsonic regions. 

The surface pressure distribution along the leeward 
plane of symmetry in the nose region is presented in Fig. 7b. 
At z = 170 in., which is the beginning of the canopy, the 
pressure increases rapidly from C p = 0.3 to » 1.0. An 
embedded subsonic pocket is formed in the canopy region 
and required three relaxation marching sweeps to develop 


x = 725" x = 925" 


x = 1200" 


Fig. 5a Shuttle Orbiter configuration and sample grids- 


vs 

1.0 


c p 

0 5 
0 


■lirO 

Fig. 5b Surface pressure distribution 
along upper symmetry plane 





Fig. 5c Canopy region pressure contours for 
x - y and z - y sectional cuts 
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the solution. The results show that the present predic- 
tion is in excellent agreement with data. Pressure contours 
on the upper symmetry plane and on the marching plane 
cross-sectional views are shown in Fig. 5c. The shock and 
expansion waves induced by the canopy can be clearly seen 
in this figure. 

Figure 5d shows the details of the orbiter geometry in 
the OMS pod region as modeled in this study. A detached 
OMS pod shock and a large subsonic pocket are formed 
in this region. Since the subsonic pocket is big and the 
Mach number is almost sero near the root of the OMS pod, 
a total of 30 relaxation marching sweeps (forward only) 
are required to give a good, converged result. Figure 5e 
presents the pressure and Mach number contours as ob- 
tained in this region. The cross-sectional pressure contours 
at x = 1080 in. and x = 1125 in. are given in Fig. 5f. The 
OMS pod shock is formed around x = 1050 in., then grows 
and finally hits the upper wing surface at x » 1095 in.. The 
chordwise pressure distributions on the upper surface of the 
wing at several span stations are compared with experimen- 
tal data in Fig. 5g. It is seen that the present calculation 
agrees with the experimental data very well over the entire 
upper surface including in the region where the OMS pod 
shock interacts with the wing surface. 



Fig. 5d Computational surface geometry of Orbiter 
at OMS pod region 


PRESSURE CONTOURS 



MACH NUMBER CONTOURS 



X 


Fig. 5e Pressure and Mach number contours 
at OMS pod region: x — y section 
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PRESSURE CONTOURS PRESSURE CONTOURS 



Fig. 5f Pressure contours at OMS pod region: 
z — y section 



Fig. 5g Orbiter upper surface pressure distribution 


?,Q —CONCLUDING JtEMARKS 

A new computational procedure has been devised to 
solve the Euler equations for three-dimensional supersonic 
inviscid flows with subsonic pockets. The method is akin 
to a simple marching procedure in portions of the flow field 
where the component of velocity normal to the local march- 
ing plane is supersonic. When this local velocity is subsonic 
(in subsonic pockets for example), a relaxation approach is 
used. The marching and relaxation strategies are both but 
variations of a unified approach to the development of fi- 
nite difference methods for this class of problems. This 
approach is based on a planar Gauss-Seidel procedure cou- 
pled with approximate factorization in the plane. Being 
an expository paper, detailed formulae are presented to aid 
the reader who would like to program the method indepen- 
dently. 


It is of interest to note the following observations: The 
method presented is not only applicable to supersonic flows 
with subsonic pockets, but is also applicable to all com- 
pressible inviscid flow regimes including entirely subsonic, 
transonic (subsonic flow with supersonic pockets), and en- 
tirely supersonic flows. By iterating in just one march- 
ing plane, a computer program based on the method pre- 
sented can be used to also solve problems that are two- 
dimensional or that can be reduced to two dimensions. 
Conical flows are examples of the latter. The same com- 
puter program can also be used to solve three-dimensional 
(and two-dimensional) unsteady flows. Thus, the unified 
approach taken is really greater in scope and applicability 
than what the title of this paper might suggest. Of course, 
the method is eminently suitable for the case of supersonic 
flows with subsonic pockets. 

The use of TVD discretizations results in a highly re- 
liable method with no artificial parameters such as coeffi- 
cients of numerical smoothing to be provided by 'the user. 
Spurious oscillations and expansion shocks are also elimi- 
nated. 

The relaxation approach used can also be used to solve 
the Navier-Stokes equations 13 . By following the unified 
methodology used to derive the present algorithm for the 
Euler equations, a unified scheme can also be derived for the 
Navier-Stokes equations. Parabolized forms of the Navier- 
Stokes equations may be solved in regions where there is 
little upstream propagating influence upon the boundary 
layer and when the flow external to the boundary layer is 
supersonic. In regions of separation, etc. where there is an 
appreciable effect of the downstream upon the upstream 
flow, and/or where the external flow is subsonic, the re- 
laxation approach may be usee). Such methods can provide 
superior replacements to current Parabolized Navier-Stokes 
solvers and can be the subject of research by us or by other 
investigators. 


15 






ACKN OW LEDGEME NT 

This study was partially supported by NASA Lang- 
ley Research Center under Contract Number NAS 1-17492 
awarded to Rockwell International Science Center. 


[12] S. R. Chakravarthy, K. Y. Szema, J. J. Gorski, 

U. C. Goldberg, and S. Osher, “Application of a New 
Class of High Accuracy TVD Schemes to the Navier- 
Stokes Equations,* AIAA Paper No. 85-0165, January 
1985. 


REFERENCES 


[1] V. Shankar, K. Y. Szema, and S. Osher, “A Conser- 
vative Type-Dependent Pull Potential Method for the 
lYeatment of Supersonic Flow with Imbedded Sub- 
sonic Region,* AIAA Paper No. 83-1887, July 1983. 

[2] K. Y. Szema and V. Shankar, “Nonlinear Computation 
of Wing- Body-Vertical Tail- Wake Flows at Low Super- 
sonic Speed,* AIAA Paper No. 84-0427, Jan. 1984. 

[3] K. Y. Szema, W. L. Riba, V. Shankar, and J. J. Gorski, 
“Full Potential Treatment of Flows Over 3-D Geome- 
tries Including Multi-body Configurations," AIAA Pa- 
per No.' 85-0272, Jan. 1985. 

[4] * S. R. Chakravarthy and D. K. Ota, “Numerical Issues 
in Computing Inviscid Supersonic Flow Over Conical 
Delta Wings," in preparation. 

[5] S. R. Chakravarthy and S. Osher, “A New Class of 
High Accuracy TVD Schemes for Hyperbolic Conser- 
vation Laws,* AIAA Paper 85-0363. 

[6] S. R. Chakravarthy and S. Osher, “Computing With 
High-Resolution Upwind Schemes for Hyperbolic Equa- 
tions*, to appear in the Proceedings of the 1983 AMS- 
SIAM Summer Seminar on Large-Scale Computations 
in Fluid Mechanics, published by American Mathemat- 
ical Society in Lectures in Applied Mathematics , Vol- 
ume 22, 1985. 

[7] P. L. Roe, “Approximate Riemann Solvers, Parameter 
Vectors, and Difference Schemes," Journal of Compu- 
tational Physics , VoL 43, 1981, pp. 357-372. 

[8] S. R. Chakravarthy, “Relaxation Methods for Unfac- 
tored Implicit Upwind Schemes," AIAA Paper 84-0165, 
1984; also to appear in the AIAA Journal. 

[9] S. R. Chakravarthy, “Numerical Implementation of 
Boundary Conditions for Hyperbolic Systems of Con- 
servation Laws,* in preparation. 

[10] S. R. Chakravarthy, “Euler Equations — Implicit 
Schemes and Boundary Conditions," AIAA Journal, 
Volume 21, Number 5, May 1983, pp. 699-706. 

[11] J. C. Townsend, D. T. Howell, I. K. Collins, and 

C. Hayes, “Surface Pressure Data on a Series of Ana- 
lytic Forebodies at Mach Numbers from 1.7 to 4.5 and 
Combined Angles of Attack and Sideslip,* NASA TM 
80062, June 1979. 



1. Report No. 2. Government Accession No. 

NASA CR-3982 


4 . Titl* l/xl Subtitle 

FULL POTENTIAL .METHODS FOR ANALYSIS/DESIGN OF COMPLEX 
AEROSPACE CONFIGURATIONS 


3. Recipient'* Catalog No. 


5. Report Oat* 

May 1986 


6. Performing Organiiation OxJe 


7. Authorial 

Vijaya Shankar, Kuo-Yen Szema, and Ellwood Bonner 


9. Performing Organization Name and Addrett 

Rockwell International Corporation. 
Los Angeles, CA 90009 


12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, DC 20546 


8. Performing Organization Report No. 


10. Work Unit No. 


11. Contract or Grant No. 

NASl-15820 


13. Type of Report and Period Covered 

Contractor Report 


14. Sponsoring Agency Code 

505-61-71-05 


15. Supplementary Note* 

Langley Technical Monitors: Noel A. Talcott, Jr., and Kenneth M. Jones 


16. Abstract 

The steady form of the full potential equation, in conservative form, is 
employed to analyze and design a wide variety of complex aerodynamic shapes. 

The nonlinear method is based on the theory of characteristic signal propagation 
coupled with novel flux biasing concepts and body-fitted mapping procedures. 

The resulting codes are vectorized for the CRAY-XMP and the VPS-32 supercomputers. 

V 

Use of the full potential nonlinear theory is demonstrated for a single- 
point supersonic wing design and a multipoint design for transonic maneuver/ 
supersonic cruise/maneuver conditions. Achievement of high aerodynamic effi- 
ciency through numerical design is verified by wind tunnel tests. Other studies 
reported here include analyses of a canard/wing/nacelle fighter geometry. 


17. Key Words (Suggested by Autbor(i)) 

Aerodynamics 
Numerical Methods 
Aircraft Design/Analys is 
Potential Theory 

18. Distribution Statement 

FEDD Distribution 
Subject Category 02 

19. Security Qaiuf. (of thri report) 

20. Security CU«if. (of thii page) 

21. No. of Pa gn 

22. 

Unclassified j 

Unclassified 


160 



Available: NASA's Industrial Applications Centers 


NASA-Langley, 1986 

























